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Numerical  Manifold  Method 


Gen-hua  Shi 

Geotechnical  Lab,  US  Army  Engineer  Waterways  Experiment  Station 
Vicksburg,  MS  39180-6199 


Abstract 

Aiming  at  global  analysis,  the  well  known  mathematical  manifold  is  perhaps  the 
most  important  subject  of  modem  mathematics.  Based  upon  mathematical  manifold,  this 
numerical  manifold  method  is  a  newly  developed  general  numerical  method.  This  method 
computes  the  movements  and  deformations  of  structures  or  materials.  The  meshes  of  the 
numerical  manifold  method  are  finite  covers.  As  the  material  domains,  the  finite  covers 
overlapped  each  other  and  covered  the  entire  material  volume.  On  each  cover,  the  manifold 
method  defines  an  independent  cover  displacement  function.  The  cover  displacement  func¬ 
tions  on  individual  covers  are  connected  together  to  form  a  global  displacement  function  on 
the  entire  material  volume. 

The  global  displacement  function  are  the  weighted  averages  of  local  independent 
cover  functions  on  the  common  part  of  several  covers.  Using  the  finite  cover  systems,  con¬ 
tinuous,  jointed  or  blocky  materials  can  be  computed  in  a  mathematically  consistent  manner. 
For  a  manifold  computation,  the  mathematical  mesh  and  physical  mesh  are  independent. 
Therefore,  the  mathematical  mesh  is  free  to  define  and  free  to  change.  As  the  mathematical 
mesh,  the  covers  can  be  moved,  can  be  split  and  can  be  easily  removed  and  added.  Mov¬ 
ing  the  covers,  the  large  deformations  and  moving  boundaries  can  be  computed  by  steps. 
By  dividing  a  cover  to  two  or  more  independent  covers  with  their  displacement  functions, 
jointed  and  blocky  materials  can  be  modeled. 

Both  the  finite  element  method  (FEM)  for  continua  and  the  discontinuous  deforma¬ 
tion  analysis  (DDA)  for  block  systems  are  special  cases  of  this  numerical  manifold  method. 
In  the  current  development  stage  of  numerical  manifold  method,  by  using  finite  cover  ap¬ 
proach,  the  extended  finite  element  method  can  compute  more  flexible  and  visible  deforma¬ 
tions  and  movements  of  joints  and  blocks. 
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Chapter  1 

General  Finite  Covers  of  Manifold  Method 

1.1  Finite  Covers  Formed  by  Mathematical  Mesh  and  Physical  Mesh 

Based  on  finite  cover  systems,  the  newly  developed  “manifold  method”  has  the  po¬ 
tential  to  meet  more  engineering  requirements.  The  term  “manifold”  here  is  a  generalization 
of  the  differential  manifold,  which  is  the  main  subject  of  differential  geometry,  topology, 
differential  topology  and  modem  algebra  of  mathematics.  The  difference  of  the  “manifold” 
here  and  the  traditional  differential  manifold  is  the  following:  the  global  functions  of  the 
differential  manifold  are  highly  differentiable  and  entirely  defined  irrelevant  with  the  cov¬ 
ers;  the  global  functions  of  manifold  here  defined  are  based  on  covers  and  only  piecewise 
differentiable,  mostly  discontinuous  on  the  contact  interfaces. 

Physically,  material  objects  often  have  different  shapes.  When  the  material  volumes 
have  fractures,  blocks  or  different  zones,  the  shapes  and  boundaries  become  more  complex. 
Under  conditions  of  large  deformation  and  moving  boundaries,  more  difficulty  occurs  be¬ 
cause  the  conventional  analytical  approximations  are  feasible  and  useful  only  in  a  local  con¬ 
tinuous  domain  which  represents  only  a  small  part  of  whole  material  volume. 

Manifolds  connect  many  individual  folded  domains  together  to  cover  the  entire  ma¬ 
terial  volume.  Then,  the  global  behavior  can  be  computed  by  functions  defined  in  local  cov¬ 
ers.  The  new  method  has  separated  and  independent  mathematical  mesh  and  physical  mesh: 
the  mathematical  mesh  defines  only  the  fine  or  rough  approximations;  as  the  real  material 
boundary,  the  physical  mesh  defines  the  integration  fields. 

The  mathematical  mesh  is  chosen  by  user,  consists  of  finite  overlapping  covers 
which  occupy  the  whole  material  volume.  Conventional  meshes  and  regions,  such  as  reg¬ 
ular  grids,  finite  element  meshes  or  convergency  regions  of  series,  can  be  transferred  to  fi¬ 
nite  covers  of  the  mathematical  mesh.  Based  on  finite  covers,  manifold  method  is  flexi¬ 
ble  enough  to  contain  and  combine  well  developed  analytical  method,  widely  used  FEM 
method  and  joint  or  block  oriented  DDA  method  in  a  unified  form. 
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FIGURE  1.1  General  covers  with  one  joint 


FIGURE  1.2  General  covers  with  two  joints 


The  physical  mesh  includes  the  boundaries  of  the  material  volume,  joints,  blocks 
and  the  interfaces  of  different  material  zones.  The  constantly  changing  water  surfaces  are 
also  part  of  the  physical  meshes.  The  physical  mesh  represents  material  conditions  which 
can  not  be  chosen  artificially. 

The  physical  cover  system  is  formed  by  both  mathematical  and  physical  meshes. 
If  the  joints  or  block  boundaries  divide  a  mathematical  cover  to  two  or  more  completely 
disconnected  domains,  those  domains  are  defined  as  physical  covers.  Therefore,  the  physi¬ 
cal  covers  are  the  subdivision  of  the  mathematical  covers  by  discontinuities.  The  manifold 
method  is  more  suitable  to  compute  large  deformations,  moving  boundaries  of  both  contin¬ 
uums  and  jointed  materials. 

In  Figure  1.1  and  1.2,  two  circles  and  one  rectangle  (indicated  by  thin  lines)  delimit 
three  mathematical  covers 

Vu  V2,  V3 

to  form  the  mathematical  mesh.  The  thick  lines  indicate  the  material  boundary  and  inner 
curved  joints.  In  Figure  1.1,  V\  is  divided  by  the  physical  mesh  into  two  physical  covers 
li,  I2,  V2  has  two  physical  covers  2i,  22  and  V3  has  two  physical  covers  3i,  32. 

Figure  1.2  shows  a  more  complex  mesh.  Mathematical  cover  V2  contains  three 
curved  lines,  but  only  two  totally  disconnected  physical  covers  2i ,  22  are  formed.  The  up¬ 
per  curve  (inside  cover  2i)  can’t  cut  through  rectangle  V2  to  form  more  physical  covers, 
therefore  cover  2i  is  a  single  physical  cover.  Similarly  since  mathematical  cover  V3  just 
intersects  the  end  of  the  upper  curves,  physical  covers  3i,  32  are  formed.  In  both  Figure 

1.1  and  1.2,  the  common  part  of  two  or  more  physical  covers  are  defined  as  “elements”  and 
marked  by  its  cover  numbers. 

Figure  1.3  shows  a  simple  but  often  useful  chain  cover  system.  This  cover  system 
is  specially  convenient  for  long  and  narrow  material  shapes. 

Figure  1.4  shows  a  DDA  block  system,  here  each  block  is  a  mathematical  cover  and 
a  physical  cover.  There  are  no  overlaps  between  any  two  covers  in  DDA  case. 

1.2  Cover  Functions  and  Weight  Functions  on  Finite  Covers 

The  normal  analytical  method  or  series  method  work  only  on  very  simple  domains 
such  as  spheres,  rectangles  and  the  functions  are  limited  to  be  highly  differentiable.  For  the 
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FIGURE  1.3  Chain  cover  system 
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manifold  method,  the  cover  displacement  functions  are  independently  defined  on  individual 
physical  covers.  Local  displacement  functions  can  be  connected  together  to  form  a  global 
displacement  function  on  the  whole  material  volume.  Since  the  joints  can  cut  one  cover  to 
more  covers  as  shown  in  Figure  1.1  and  Figure  1.2,  the  functions  are  disconnected  on  the 
two  sides  of  these  joints.  The  global  displacement  function  is  general  and  flexible  enough  to 
represent  the  wide  variety  of  continuous  or  discontinuous  materials  located  within  moving 
boundaries. 

The  cover  functions  u,(x,  y)  defined  on  physical  cover  Ui 


ui(x,y)  (*£? y)  €  Ui  (!•!) 

can  be  constant,  linear,  high  order  polynomials  or  locally  defined  series.  These  cover 
functions  are  connected  together  by  the  weight  functions  u>j(x,  y ) 


Wi(x,y)  >0  ( x,y )  €  Ui 

Wi(x,y)  =0  (x,y)£Ui  (1.2) 

with 

X]  Wj(x,y)  =  l.  (1.3) 

(x,y)€Uj 

The  meaning  of  the  weight  functions  Wi(x,  y)  is  weighted  average,  which  is  to  take 
a  percentage  from  each  cover  function  Ui(x,  y)  for  all  physical  covers  Ui  containing  (x,  y) 

Using  the  weight  functions  Wi (x,  y)  a  global  function  F(x,  y)  on  the  whole  physical 
cover  system  is  defined  from  the  cover  functions 


n 

F(x,y)  =  ^ttfj(x,i/)uj(x,y).  (1.4) 

1=1 

Figure  1.5  is  a  one  dimensional  example,  there  are  three  physical  covers 


Ui  =  AiA2,  U2  =  BXB 2,  U3  =  CXC2 


l 


cover  U2 


FIGURE  1.5  1-d  general  covers 
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Ul(x)  = 

^■3^4  » 

x  e  Ui 

U2{x)  = 

B3B4  , 

x  €  U2 

u3(x)  = 

C3C4, 

x  e  Uz 

toi(a;)ui(x)  = 

ri.3ri.5A2, 

x  e  U\ 

w2(x)u2(x)  = 

BiB^BeB2, 

x  €  U2 

w3(x)u3(x)  = 

C2C5C4, 

xeUz 

The  global  function  F(x )  is 

n 

F(x)  =  ^2  wi(x)ui(x)  =  AiA^B^B^C^Ci. 

i=  1 


1.3  Global  Functions  on  Continuous  and  Discontinuous  Materials 

For  material  analysis,  four  basically  different  methods  are  often  used.  In  the  order  of 
their  development,  analytical  solutions  are  the  earliest,  then  came  finite  difference,  the  finite 
element  method,  and  most  recently  the  distinct  element  method  and  the  discontinuous  de¬ 
formation  analysis.  The  analytic  approach  has  simple  material  boundary  such  as  rectangles 
or  circles  and  has  no  meshes.  The  finite  difference  uses  grid  meshes  with  equal  spaces  and 
as  such,  is  more  general  than  the  analytic  solution  method.  The  finite  element  method  was  a 
revolution,  it  shifts  from  differential  equations  to  integral  equations,  from  the  smooth  func¬ 
tions  to  the  piece  wise  smooth  functions.  The  meshes  of  finite  element  method  can  give 
results  of  generally  shaped  continuous  materials.  The  latest  distinct  element  method  and 
discontinuous  deformation  analysis  method  are  for  block  systems  which  are  totally  discon¬ 
tinuous.  The  displacement  functions  of  distinct  element  method  and  discontinuous  defor¬ 
mation  analysis  are  defined  for  individual  blocks  of  general  shape  which  are  completely 
disconnected  from  block  to  block. 

As  a  one  dimensional  example.  Figure  1.6  shows  the  accuracy  of  different  methods 
when  approaching  a  natural  function  (thin  curves)  which  is  discontinuous  at  a  point.  The 
thick  smooth  curve  of  Figure  1.6  a)  is  the  approximation  from  the  analytical  and  finite  differ¬ 
ence  methods.  The  thick  piecewise  smooth  segments  of  Figure  1.6  b)  are  the  approximation 
from  FEM.  The  one  dimensional  elements  are 


Xq-Ti,  X\X2,  X2X3,  X3X4,  x±x5. 
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The  disconnected  segments  of  Figure  1.6  c)  are  the  approximations  from  DEM  and 
DDA  methods.  The  one  dimensional  blocks  are 

yoZl,  2/1^2,  J/2®3,y3Z4,  S/4Z5- 

which  have  more  unknowns  than  the  previous  methods. 

Figure  1.6  d)  and  Figure  1.7  show  the  approximation  of  the  manifold  method.  There 
are  seven  one  dimensional  physical  covers 


Ui=x0xi,  U2  =  x0x2,  U3  =  xix3,  U4  =  x2x3, 

Ub=yzx4,  Ue  =  y3x5,  U7  =  x4x5 

Since  the  natural  function  has  a  jump  at  the  point  x3  =  y3,  the  mathematical  cover 
x2x4  was  split  to  two  physical  covers  U4  =  x2x3  and  f/5  =  y3x4. 


mi(.r)iti(.'c)  = 

AqX  1, 

x  e  U\ 

tw2(ar)tx2(a:)  = 

XqA^X2, 

x  e  u2 

w3(x)u3(x)  = 

X\  A2x3, 

x  <E  Uz 

tU4(x)tt4(a;)  = 

x2A3 , 

x  £  U4 

w5(x)u5(x)  = 

B3x4, 

x  €  CT 5 

w6{x)u6(x)  = 

yj,A4xh , 

x  6  U6 

iu7(x)u7(x)  - 

x4A$, 

x  G  U7 

The  global  function 

7 

F(x)  =  Wj(x)uj(x)  =  Aq A\  A2A3B3A4  A5  (1-5) 

j=i 

is  very  close  to  the  original  natural  curve.  The  global  displacement  functions  of  the 
manifold  method  are  capable  of  representing  large  deformations  of  fractured  or  blocky  ma¬ 
terials  until  the  ultimate  damage  stage  in  a  unified  mathematical  form. 

For  the  discontinuous  deformation  analysis,  the  material  body  is  simply  individual 
blocks.  Each  block  is  a  mathematical  cover,  and  each  mathematical  cover  is  a  physical 
cover.  The  mathematical  mesh  and  physical  mesh  are  the  same  where  all  covers  are  not 
overlapped.  Therefore  the  discontinuous  deformation  analysis  is  the  totally  discontinuous 
case  of  manifold  method. 


11 


1 .4  Elements  As  The  Common  Part  of  Overlapped  Covers 

For  two  dimensional  manifold  computation,  the  cover  displacement  functions 
Ui(x,y)  and  Vi(x,  y )  are  defined  on  cover  The  global  displacement  functions  u(x,  y) 
and  v(x,  y)  on  the  whole  material  body  are  two  global  functions: 

The  cover  displacement  functions  ui(x,  y)  and  Vi(x,y)  defined  on  physical  cover 

Ui 


Ui(x,y)  {x,y)  e  Ui 
Vi(x,y )  ( x,y)eUi 

can  be  constant,  linear,  high  order  polynomials  or  locally  defined  series.  These  cover 
displacement  functions  are  connected  together  by  the  weight  functions  u>i(x,  y) 


Wi(x,y)>  0  ( x,y)6Ui 

Wi(x,y)  =0  (x,y)£Ui  (1.6) 


with 

Y  Wj{x,y)  =  l 

( x,y)eUj 

Using  the  weight  functions  Wi{x,y )  a  global  displacement  function  u(x,y)  and 
v(x ,  y)  on  the  whole  physical  cover  system  is  defined  from  the  cover  displacement  func¬ 
tions 


n 

u(x,y)  =  Ywi(x^y)udx^y) 

l  —  l 
n 

■u(-uy)  =  YWi(x,y}Vt(x'y} 

i  —  1 


Each  common  domain  of  physical  covers  is  defined  as  a  element  e.  From  formula 
(1.6),  in  each  element,  weight  function  wt(x,  y)  has  a  analytical  representation,  which  is 
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either  constant  0  or  a  differentiable  elementary  function.  Therefore,  global  displacement 
function  (u,  v )  has  analytical  representation  in  each  element. 


The  cover  displacement  functions  can  be  given  in  the  following  form: 


(  ^yA 

\v{x,y)  ) 

=  PHx,y)]  [Di] 


( tn{x,y)  t12(x,y)  tn(x,y ) 
V2i(z,y)  t22(x,y)  t23(x,y) 


tlm-l(x,y) 

^2m-l(z,y) 


tim(x,y)\ 
hm{x,y)  ) 


d il  \ 

di2 

di  3 

di4 


I  dim-1  I 

V  dim  / 

(1.7) 


where  the  subscript  “i”  represents  the  i-th  cover. 

If  a  element  e  is  not  in  cover  Ui,  T,(x,  y)  =  0. 

Ti{x,y)  ^ 0  e  <E  Ut 
Ti(x,y)  =  0  e  g  Ui 


1.5  Possible  Displacement  Functions  for  General  Covers 

On  cover  Ui,  ( Ui(x ,  y),  V{(x,  y))  is  the  displacements  of  point  (a:,  y)  on  x  any  y  di¬ 
rection  respectively,  The  cover  displacement  functions  (ui(x,  y),Vi(x,  y))  often  take  one  of 
the  following  forms: 

the  constant  function  on  cover  Ui 

Ui(x,y)  =  dn 

Vi{x,y)=dt2  (1.8) 

the  complete  first  order  approximation  on  cover  Ut 

Ui(x,y)  =  dn  +  di3x  +  dihy 

Vi(x,y)  =  di2  +  dux  +  di6y  (1.9) 
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the  complete  second  order  approximation  on  cover  Ui 


Ui(x,y)  =  dn  +  di3x  +  di5y  +  di7x 2  +  di9xy  +  dmy2 

Vi{x ,  y)  =  di2  +  dux  +  di6y  +  disx 2  +  dil0xy  +  du2y2  (1.10) 


or  the  general  series  form  U{(x,y )  and  vl(x}  y)  on  cover  Ut 


f  Ui(x,y)  1  _ 
1  uf(a,y)  /  ' 


( fi(x,y)  0 

V  0  /i(*>y) 


f2{x,iy)  0 

0  h(x,v) 


( 


fm{x,y)  0  \ 

0  fm(x,y) ) 


du  N 
di2 
du 
du 


I  di2m  —  \  I 
\  di2m  / 
(1.11) 


1.6  Coefficient  Matrix  of  Equilibrium  Equations  for  General  Covers 

Assume  the  number  of  physical  covers  is  n,  and  there  are  2m  unknowns  in  each 
physical  cover. 


/ 


Di  = 


du 
di  2 
diz 
du 


\ 


i  =  1,2, , . . ,  n  . 


I  di2rn  —  l  I 
\  d{2m  / 


the  total  potential  energy  has  the  form 
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/An 

K 12 

co 

•  •  •  Fin  ^ 

(Dl\ 

K-n 

I<22 

K2Z 

...  K 2n 

d2 

Dj 

DJ 

...  DTn) 

K3 1 

Kz2 

Kzz 

...  I<Zn 

Dz 

\I<nl 

I<n2 

I<nZ 

•••  I<nn) 

\Dn) 

+  (DJ  Dj  Dj  ...  Dl) 


(F'\ 

F2 

Fz 

\Fn/ 


+  C 


(1.12) 


Because  each  cover  has  2m  degrees  of  freedom,  each  submatrix  I\l]  in  the  coefficient  matrix 
given  by  equation  (1.12)  is  a  2m  x  2m  matrix.  D,  and  Fi  are  2m  x  1  submatrices,  where 
Di  represents  the  displacement  variables  (dn  dt2  d{3  . . .  di2m)T  of  physical  cover  i. 

From  the  formulation  of  II,  the  formula  (1.12)  can  be  written  as  a  symmetric  representation, 


Ki 


The  equilibrium  equations  are  derived  from  minimizing  the  total  potential  energy 
II.  The  i-th  row  of  following  equation  (1.13)  consists  of  2m  linear  equations 

5X1 

——  =0,  r  =  1, 2, 3,4, ... ,  2m, 
ucl{r 

where  dir  is  the  displacement  variable  of  cover  i.  The  matrix  of  obtained  simultaneous  equi¬ 
librium  equations  is  same  as  the  matrix  of  quadratic  form  (1.12): 


/An 

K 12 

A' 13 

•  K\n  \ 

(Dl\ 

fFl\ 

A21 

A'22 

I< 23  • 

■  I<2n 

d2 

f2 

Kzi 

A'32 

Kzz  • 

■  Kzn 

Dz 

— 

Fz 

\Kn  1 

Kn  2 

Knz 

•  KnJ 

\Dn) 

\Fn) 

For  material  analysis,  Ft  is  the  loading  on  cover  i  distributed  to  the  2m  displacement  vari¬ 
ables.  Submatrices  [//,,]  depend  on  the  material  properties  of  cover  i  and  [A',y] ,  where  i  /  j 
is  defined  by  the  overlapping  or  contact  between  cover  i  and  cover  j. 
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Chapter  2 

Element  Matrices  of  General  Finite  Covers 


2.1  Stiffness  Matrix  for  General  Covers 


For  the  manifold  method  with  general  covers,  the  elements  are  the  common  domains 
intersected  by  the  physical  covers.  The  integration  domains  of  the  stiffness  matrices  are  the 
whole  elements  which  are  the  intersection  of  several  physical  covers. 


Same  as  FEM  method,  the  relationship  between  stress  and  strain,  is  given  by 


where 


and  E,  u  are  Young’s  modulus  and  Poisson’s  ratio  respectively.  And 


(2.1) 


u(x,y) 

v(x,y) 


1=1 


=  ([Tl]  [T2]  ...  [: Tn ]) 

\ 


/  du(x,y) 

I  dx 

dv(x,y) 
dy 

du(x,y)  dv(x,y) 
dy 


dx 


{Di} 

V{S}/ 


(2.2) 


[Ti{x,y)\  = 


ftn(x,y)  t12(x,y)  tn(x,y)  ...  t  lm- l(^?y)  tlm  (x,y)\ 

\hi{x,y)  t22(x,y)  t23(x,y)  ...  t2m-i(x,y)  t2m(x,y) ) 
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Then 
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=  \  jj  {D}T\B)T\E\[B]{D}dxdy 

=  \{d)t  [fjmm** 

=  V}r(s'[B]T[£][Bl){B}, 


{£>} 


where  Se  is  the  area  of  that  element. 


Therefore, 


Se[B}T[E}[B]  =  Se 


/*£\ 

mT 

[Bz)T 


'  [Bn]3' 


is  the  element  stiffness  matrix. 


(2.7) 


[B]  ( [Bi]  [B2|  [B,]  ...  [B„]),  (2,8) 


Then 


Se[Br)T{E}[B9}  -4  [Krs],  r,  s  =  1,2,3, ...  ,n. 


2.2  Initial  Stress  Matrix  for  General  Covers 


Following  the  time  sequence,  the  manifold  method  computes  step  by  step.  The  com¬ 
puted  stresses  of  previous  time  step  will  be  the  initial  stresses  of  the  next  time  step.  There¬ 
fore  the  initial  stresses  are  essential  for  manifold  computation. 


o 

X 


For  the  element  e,  the  potential  energy  of  the  initial  constant  stresses  {<x°}  = 

^0  _0  \T  • 

°y  Txy  )  1S 


n„ 


+  eya°y  +  7rj,r°y)  dx  dy 


(  a°x\ 

ev  Ixy )  oj  )  dx  dy 

\  '  / 

\  xy  / 


JJ  {D}T[B}T{<r°e}dxdy 

S‘{D}T[B]T{a°e}, 


(2.9) 
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where  Se  is  the  area  of  that  element.  Therefore, 


is  the  load  matrix. 


-Se[Be]T{a°e}  =  -Se 


\{B'n]TJ 


(2.10) 


Then 


-Se[Br)T 


{-Fr})  r  —  1,2,3, ...  ,n. 


2.3  Point  Loading  Matrix  for  General  Covers 

Different  from  ordinary  FEM  method,  a  load  point  can  be  any  point  in  its  element. 

rp 

The  point  loading  force  (  Fx  Fy)  acts  on  point  (x,  y)  of  element  e.  And  the  displace¬ 
ments  on  force  point  (x,  y)  is 


v-/ 

The  potential  energy  due  to  the  point  loading  is 


n„  =  -(Fxu(x,  y )  +  Fyv(x,  y)) 

=  -(u(x,y)  v(x,y))^jjf^ 

=  -{c}T[r(x,j/)]T(^. 

Therefore, 

/Pi]T\ 

[ t2)t 

mT 

\[Tn]T) 


(2.11) 


(2.12) 
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is  the  load  matrix. 


Then 


{i^r  },  r  —  1)  2, 3, . . . ,  n. 


2.4  Body  Loading  Matrix  for  General  Covers 

T 

Assuming  that  (fx  fy)  is  the  constant  body  force  loading  acting  on  the  material 
area  of  element  e. 

The  potential  energy  due  to  the  body  loading  is 


n,„  = 


JJ^(fxu(x,  y)  +  fyv(x,  y ))  dx  dy 
jj r  ( «(*,  y)  v(x,  dx  dy 


-{D}t  JjjTfdxdy 


fx 

fy 


Therefore, 


\s !at][T  dxdy 

(«- 

IL 

/[T1(x,y)]T\ 
[T2(x,y)}T  1 
[T3(x,y)]T 

dx  dy 

'  [Tn(x,  y)]T  ) 

Jh[Ti{x,y)]T  dx  dy\ 
iL[T2(a:,y)]T^  dy 
jfA[Ux,y)}Tdxdy 


\fJA{Tn(x,y)}T  dx  dy) 


IJATi(x^y)]T  dxdy  \ 
jfA  {T2{x,y)}T  dxdy 
JfA[T3(x,y)}T  dx  dy 

\)JA[Tn{xyy))T  dx  dy) 


(2.13) 


(2.14) 
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is  the  body  loading  matrix.  And 


Then 


JJ  [Tr(x,y)f 


dx  dy 


{■^r})  r  —  1, 2, 3, . . . ,  n. 


2.5  Inertia  Force  Matrix  for  General  Covers 

Inertia  force  matrix  is  equivalent  to  mass  matrix  of  FEM.  This  matrix  is  the  most 
important  matrix  of  manifold  method.  Giving  a  small  time  steps,  inertia  force  matrix  will 
control  the  movements  and  the  stability  of  all  points  of  the  whole  material  volume.  In  each 
time  step,  the  displacements  should  be  small  enough  to  have  the  final  result  independent 
from  the  choosing  of  the  specific  time  steps. 

Considering  the  current  time  step,  denote  ( u(x,  y,  t)  v(x,  y,  t) )  as  the  time  de¬ 
pendent  displacements  of  any  point  (x,  y)  of  element  e  and  M  as  the  mass  per  unit  area. 
The  force  of  inertia  per  unit  area  is 


/  fx(x,y,t ) 


(  u(x,y,t ) 
\v(x,y,t) 


)  - 


The  potential  energy  of  the  inertia  force  of  element  e  is 


(2.15) 


Hr  =  -  JJ  ( u(x,  y,  t)  v{x,  y,t))(^  ^  dx  dy 


=  JJ^M{u(x,y,t)  v(x,  y,  t ) )  [T]  — —  dx  dy. 


(2.16) 


Assume  {T>(0)}  =  {0}  as  the  element  displacements  at  the  beginning  of  the  time 
step,  {D(  A)}  =  {D}  as  the  displacements  at  the  end  of  the  time  step,  A  as  the  time  interval 
of  this  time  step.  Then 


{D}  =  {D(A)}  =  {D(0)}  +  A 


d{D(0)}  A2  <92{£>(0)} 


dt 


+ 


dt2 


Ad{pm  |  A2a2{D(o)} 


dt 


dt 2 


(2.17) 


22 


=  —{D}~  1^221  =  A{J9}  _  l{y(°)}  (2.18) 

dt 2  A2  X  A  dt  A2X  5  Ax  v  v  ’ 


where 


mo)}  = 


dim) 

dt  ’ 


(2.19) 


is  the  velocity  at  the  beginning  of  the  time  step.  Then  the  potential  energy  becomes 


n,  =  M{D( A)}t  \^j j[T(x,y)]T\T(x,y)\dxdV 


(A{B(A)}  -  |{V(0»)  ■ 


(2.20) 


Therefore, 


-\JJjT(x,y)}T{T(x,y)]dXJy\  (^(»(A))-*{F(0))).  (2.21) 


is  the  load  matrix.  Thus, 


2  M 
A2 


\JJjr^x'  y)}T\-T(x'  y )1 dx  dy 


r 

/[Ti)T\ 

- 

2M 

A2 

JL 

mf 

py1' 

aril  in\  mi  . 

.  [Tn] )  dx  dy 

V[T„fJ 

- 

[*], 


(2.22) 


and 


2A4 

A 

2M 
''  A 


[JJjT(x ,  y)]r[T(x,  y)]  da:  dy 


WO)} 


“ 

JL 

inf 

pyT 

(Pi]  Pi]  [T3]  . 

.  [T„] )  dx  dy 

[TJT  / 

WO)}  ->  {F}. 


(2.23) 
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Then 


2  M 

A2 


IL 


[Tr(x,y)]T[Ta(x,y)]  dxdy 


V)  s  1,  2,3, . . .  ,n. 


[Krs],  (2.24) 


2  M 
A 


[Tr{x,y)]T[T,(x,y)}dx  dy 


{V»(0)} 


r,  5  =  1,2,3, . . . ,  n. 


where 


W°)}  =  * 


/'".(OA 

U(° ))' 


(2.25) 


2.6  Fixed  Point  Matrix  for  General  Covers 

As  a  boundary  condition,  some  of  the  elements  are  fixed  at  specific  points.  The  con¬ 
straint  can  be  applied  by  using  two  very  stiff  springs.  Assume  the  fixed  point  is  (x,  y)  at 
element  e  with  the  displacements 


( u(x,y)\  f0\ 

\v(x,y)J  V0/ 

There  are  two  springs  which  are  along  the  x  and  y  directions  respectively.  The  stiff¬ 
ness  of  the  springs  is  p.  The  spring  forces  are 

f  f A  _  f~pu{x,y)\ 

\fyj  \-pv{x,y)J' 


The  strain  energy  of  the  spring  is  II/,  then 


n/  =  7,{u(x,y)2  +v(x,y)2) 

=  !< «(*•»)  «<*•»»  (“{*;;} 


(2.26) 
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Then 


p[T(x,y))T\T(x,y)] 


=P 


/[Ti(*,y)}T\ 

PM*,y)]T 

[T3(x,y)]r 


([Ti(xty)\ 


Pifoy)]  [T3(x,y)} 


[■ Tn(x,y )] ) . 


\{Tn(x,y)]T/ 


(2.27) 


is  the  [K]  matrix.  Therefore, 


p[Tr(x,y)]T[Ts(x,y)] 


[Krs\,  r,s  =  1,2,3,...,  n. 
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Chapter  3 

Entrance  Theory  on  Contacts  for  General  Covers 


3.1  Definition  of  Contacts  for  General  Covers 

Thus  far,  only  individual  covers  and  elements  were  considered.  It  is  necessary  to 
connect  the  individual  discontinuous  boundaries  into  a  system.  For  the  movements  of  dis¬ 
continuous  boundaries,  no  tension  and  no  penetration  must  be  satisfied  between  two  contact 
sides.  The  following  theory  of  kinematics  for  material  is  based  on  contacts.  Finding  the  con¬ 
tacts  in  each  time  step,  allying  stiff  springs  on  contacts,  the  discontinuous  displacements  can 
be  computed. 

The  time  steps  can  be  chosen  small  enough  so  that  the  displacements  of  all  points 
in  the  whole  material  body  are  less  than  a  pre-defined  limit  p.  Therefore,  the  displacements 
can  be  small  enough  so  that  the  displacements  (u(x,  y),  v{x,  y )),  the  rotations  r(x,  y),  the 
deformations  (ex  ey  ^xy)  can  be  accurately  represented  as  linear  functions  of  the  cover  un¬ 
knowns  [D{] 


(u(x,  y),v(x,y))  =  ]P[Tf(x,y)]{A} 
i=  1 

n 

rix,y)  =  ^[i?,(a:,y)]{A} 

t=  1 

/  €x  \  n 

(  £y  )  =  ^[5,(x,y)]{Di}  (3.1) 

\  7 xy  /  l=l 


Based  on  the  small  step  displacements,  the  contacts  are  defined  in  the  beginning  of 
each  time  step.  Each  contacts  are  formed  with  two  sides.  All  of  the  pair  of  two  sides  that  are 
possible  to  contact,  penetrate  or  entrance  from  one  to  another  at  the  end  of  the  time  step  are 
defined  as  contacts.  Since  practically  there  are  no  penetrations  on  the  two  sides  that  allowed 
on  the  contacts,  so  the  contacts  are  merely  the  entrance  positions. 
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There  are  two  kinds  of  contacts:  angle  to  edge  and  angle  to  angle.  The  edge  to  edge 
contacts  can  be  transferred  to  angle  to  edge  contact.  (See  Figure  3.1)  Assume  the  maximum 
step  edge  rotation  is  8,  there  are  the  criteria  for  the  contacts: 

[1]  for  angle  to  edge  contact,  the  minimum  distance  of  the  angle  vertex  to  the  edge  of 
the  contacts  is  less  than  2 p; 

[2]  for  angle  to  angle  contact,  the  minimum  distance  of  two  angle  vertices  of  the  contacts 
is  less  than  2 p;  (See  Figure  3.1) 

[3]  for  angle  to  edge  contact,  when  the  angle  vertex  translates  to  the  edge  without  rota¬ 
tion,  the  maximum  overlapping  angle  of  the  angle  and  the  edge  is  less  than  28; 

[4]  for  angle  to  angle  contact,  when  the  angle  vertex  translates  to  the  vertex  of  other 
angle  without  rotation,  the  maximum  overlapping  angle  of  the  two  angles  is  less  than 
28;  (See  Figure  3.2) 

The  computer  program  sets  the  maximum  step  edge  rotation  8  =  1.5  degrees. 

Figure  3.2  shows  two  complex  blocks  with  many  edges  and  angles,  under  this  con¬ 
tact  criteria,  there  are  only  two  contacts  even  if  the  distance  limit  2 p  is  lager  than  block 
diameters. 

The  common  sense  requirements  of  no  penetration  and  no  tension  indeed  are  in¬ 
equalities.  Coloumb’s  friction  law  and  limited  tension  are  also  inequalities.  From  formula 
(3.1),  the  inequalities  can  be  simplified  to  linear  inequalities  of  There  are  three  kinds 
of  linear  inequalities: 

[1]  no  penetrations  in  contacts, 

[2]  tension  force  less  than  tension  strength  in  contacts, 

[3]  Coloumb’s  friction  law  in  contacts. 

In  the  computation,  an  angle  to  angle  contact  will  be  transferred  to  one  or  two  angle 
to  edge  contacts.  Based  on  the  orientation,  the  following  formula  is  used  to  judge  the  angle 
to  edge  contact  penetration. 
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Assume  Pi  is  a  point  before  deformation  which  moves  to  point  P{  after  deforma¬ 
tion;  P2P3  is  the  entrance  line  and  (xj,  yt)  and  (u,-,  vt)  are  the  coordinates  and  displacement 
increments  of  P,  ,  i  =  1, 2, 3  respectively.  If  points  Pi ,  P2  and  P3  rotate  in  the  same  sense 
as  the  rotation  of  ox  to  oy  (see  Figure  3.3  and  Figure  3.4),  then  P{  has  passed  line  P2P3  and 
is  stated  by  the  inequality: 


A  = 


1  Xi  +  U\ 

1  x2  +  u2 

1  x3  +  lt3 


2/i  +  vi 
2/2 

2/3  +  ^3 


<  0. 


(3.2) 


This  simple  formula  is  still  correct  even  when  these  three  points  move  simultane¬ 
ously. 


3.2  Entrance  Lines  of  Contacts  for  General  Covers 

The  kinematics  can  be  imposed  on  the  global  equations  by  adding  stiff  springs  to 
lock  the  movement  in  one  or  two  directions.  The  theory  of  block  and  joint  kinematics  will 
decide  where  and  when  to  put  the  stiff  springs  so  that  the  movements  of  blocks  and  joints 
are  basically  the  same  as  in  the  real  space.  For  the  kinematics,  the  formulae  of  penetration 
judge  and  penetration  lock  have  to  be  consistent  to  ensue  the  convergency  of  open-close 
iterations. 

For  the  contacts,  the  “entrance”  lines  can  be  defined.  A  contact  between  two  convex 
angles  is  shown  by  Figure  3.5,  where  the  two  thick  lines  are  the  entrance  lines.  Penetration 
will  occur  if  the  two  entrance  lines  are  passed  by  the  vertices  of  the  other  angles  simultane¬ 
ously. 

[1]  For  the  angle  to  angle  contact,  if  both  angles  are  less  than  180°,  the  two  entrance 
lines  are  defined  according  to  the  following  table: 


two  entrance  lines 

a  <  180° 

P  <  180° 

OEz 

oe2 

a  <  180° 

P  >  180° 

oe3 

oea 

a  >  180° 

P  <  180° 

OE\ 

oe2 

a  >  180° 

P  >  180° 

OE\ 

oea 

[2]  For  the  angle  to  angle  contact,  if  a  angle  is  larger  than  180°,  the  two  entrance  lines 
are  the  two  edges  of  the  angle  greater  than  180°. 
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FIGURE  3.6  Entrance  lines  of  parallel  edges 


[3]  For  the  angle  to  edge  contact,  the  only  entrance  line  is  the  edge. 

[4]  For  the  edge  to  edge  contact,  the  two  edges  are  parallel,  there  is  only  one  entrance 
line  in  this  case,  the  entrance  line  can  be  any  of  the  two  edges. 

As  showing  by  Figure  3.6,  the  definition  of  entrance  line  is  still  correct  even  if  edges 
of  two  sides  of  the  contact  are  parallel  or  slightly  penetrated. 

The  penetration  judgment  in  different  contacts  are  as  following: 

[  1  ]  For  an  angle  to  angle  contact,  if  both  angles  are  less  than  180° ,  the  two  entrance  line 
have  to  be  passed  simultaneously  by  the  vertex  of  another  angle. 

[2]  For  an  angle  to  angle  contact,  if  a  angle  is  larger  than  180°,  one  of  the  two  entrance 
lines  has  to  be  passed  by  the  vertex  of  another  angle. 

[3]  For  an  angle  to  edge  contact,  the  only  entrance  line  bas  to  be  passed  by  the  angle 
vertex. 

For  an  angle  to  angle  contact,  penetration  happens  when  there  are  no  vertex  entering 
the  opposite  angle  as  shown  in  Figure  3.7.  Therefore  the  fact  that  if  a  vertex  is  in  other  angle 
is  not  a  criteria  of  penetration.  The  penetration  is  related  to  the  entrance  lines. 


FIGURE  3.7  Using  of  entrance  lines 
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In  order  to  prevent  the  penetration  of  the  two  sides  of  the  contacts,  stiff  springs  are 
applied  to  some  entrance  lines: 

[1]  For  angle  to  angle  contact,  if  both  angles  are  less  than  180°,  a  stiff  spring  is  attached 
between  the  vertex  and  its  entrance  line  which  has  been  passed  first  or  has  smaller 
passing  distance  d.  For  this  type  of  contact,  only  one  stiff  spring  is  applied. 

[2]  For  angle  to  angle  contact,  with  a  angle  larger  than  180°,  if  any  one  of  the  two  en¬ 
trance  lines  has  been  passed  by  the  corresponding  vertex  of  another  angle,  a  stiff 
spring  is  applied  along  the  normal  to  this  edge.  If  two  entrance  lines  have  been 
passed,  two  stiff  springs  are  applied  to  the  two  entrance  edges. 

[3]  For  angle  to  edge  contact,  one  stiff  spring  is  applied  to  the  entrance  edge. 

The  entrance  distance  d  from  Pi  to  P2  P3  can  be  computed  by  the  following  formula, 

1  xi  +ui  yi  +  vi 

1  X2  +  U2  2/2  +  v2  <  0.  (3-3) 

1  Xz  +Uz  2/3  +  Vz 


l  = 

If  a  contact  is  open  in  the  beginning  of  a  time  step  and  is  closed  in  the  end  of  the 
same  time  step,  the  entrance  time  and  position  in  this  time  step  can  be  computed.  Assume 
t  =  0  in  the  beginning  of  the  time  step,  and  t  =  1  in  the  end  of  the  time  step,  and 


1  X\  T  tu i  2/1  +  tvi 

A  (t)  =  1  X2  +  tu2  2/2  +  tv 2 

1  xz  T  tuz  2/3  +  tv3 

then 

A(0)  >  0  A(l)  <  0 

The  entrance  time  to  satisfies  equation 

1  x i  T  toll!  2/1  +  toVi 

A(fo)  =  1  x2  +  fou2  2/2  +  tov2  =  0 

1  Xz  T  t0uz  2/3  +  *or>3 
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(3.4) 


Neglecting  the  second  order  infinite  small,  t0  can  be  computed  by  simpler  formula: 


to 


1  yi 

1  x2  y2 
l  yz 


The  contact  position  (xo,  yo)  is 


1 

Ui 

2/1 

1 

U2 

yz 

1 

U3 

yz 

l 

Vi 

l 

X2 

V2 

l 

xz 

vz 

(  xo  \  _  (  x i  +  “i*o  \ 

\Vo  )  \Vi  T “l *o  J 


(3.5) 


(3.6) 


These  formulae  are  still  correct  even  when  the  three  points  move  simultaneously. 


3.3  Contact  Transfer  to  Next  Time  Step  for  General  Covers 

The  computation  of  manifold  method  follows  time  steps.  The  closed  contact  points 
should  go  to  next  time  step  and  find  new  representing  contacts.  The  closed  contacts  of  the 
previous  time  step  will  be  transferred  to  the  next  time  step,  if  the  contacts  are  found  in  the 
same  contact  position.  The  entrance  lines  will  be  transferred  in  case  it  is  possible.  The  angle 
to  angle  contacts  and  angle  to  edge  contacts  have  different  contact  parameters  and  different 
stiff  springs.  The  same  contact  may  also  has  different  contact  parameters  in  different  contact 
state. 


The  springs  of  closed  angle  to  angle  contacts: 

[1]  normal  stiff  spring, 

The  springs  of  closed  angle  to  edge  contacts  in  the  sliding  mode: 

[1]  normal  stiff  spring, 

[2]  pair  of  shear  sliding  forces, 

The  springs  of  the  closed  angle  to  edge  contacts  locked  in  both  normal  and  shear 
directions 

[1]  normal  stiff  spring, 

[2]  shear  stiff  spring, 
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The  closed  angle  to  angle  contacts  have  the  following  contact,  parameters: 

[1]  normal  forces, 

[2]  normal  displacements. 

The  closed  angle  to  edge  contacts  in  the  sliding  mode  have  the  following  contact 
parameters: 

[1]  normal  forces, 

[2]  normal  displacements, 

[3]  one  of  the  two  possible  sliding  directions, 

[4]  contact  positions  on  the  edges, 

The  closed  angle  to  edge  contacts  locked  in  both  normal  and  shear  directions  have 
the  following  contact  parameters: 

[  1  ]  normal  forces, 

[2]  normal  displacements, 

[3]  shear  forces, 

[4]  shear  displacements, 

[5]  contact  positions  on  the  edges. 

Transferring  the  contacts  to  next  step,  an  angle  to  angle  contact  may  transferred  to 
an  angle  to  edge  contact.  Also  an  angle  to  edge  contact  may  transferred  to  angle  to  angle 
contact.  As  transferring  to  different  kind  of  contacts,  some  contact  parameters  may  not  be 
needed. 

3.4  Normal  Contact  Matrices  for  General  Covers 

Assume  Pi  is  a  vertex,  P2P3  is  the  entrance  line  and  (a t/*)  and  («*,  Vk)  are  the 
coordinates  and  displacement  of  P&,  k  —  1, 2, 3  respectively.  If  points  Pi ,  P2  and  P3  rotate 
in  the  same  sense  as  the  rotation  of  ox  to  oy  (see  Figure  3.3),  then  the  distance  d  from  Pi  to 
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line  P2P3  is 


A  1  1  xi+ui  Vi+vi 

rf=y  =  j  1  x2+u2  y2  +  v2  ,  (3.7) 

1  X3  +u3  J/3  +  V3 

l  = 

d  should  be  negative  if  Pi  passed  edge  P2  P3 
The  step  displacements 

(■ Ui,Vi )  i  =  1,2,3 

is  small  as  the  results  of  small  time  step.  The  contact  distance  y  is  small  from  the  definition 
of  the  contacts.  Then 

1  ®i  iji  1  u\  yi  1  xi  Vi 

So  —  1  x2  V2  ,  1  u2  3/2  ,  1X2  V2 

1  X3  2/3  1  u3  y3  1  2:3  v3 

are  the  first  order  infinite  small,  and 

1  U\  V\ 

1  u2  v2 
1  u3  v3 

is  the  second  order  infinite  small. 


1  Xi  +  ux  yi+  vi 

1  X2+U2  y2  T  v2 
1  x3  +u3  J/3  +  v3 


1 

Xi 

m 

1 

U\ 

y\ 

1 

Vl 

1 

Ui 

Vl 

1 

J/2  + 

1 

U2 

V2  + 

1 

*2 

v2  + 

1 

u2 

V2 

1 

Z3 

VZ 

1 

^3 

yz 

1 

^3 

^3 

1 

^3 

V3 

Neglecting  the  second  order  infinite  small,  there  is 


1 

Ui 

yi 

1 

Xi 

Vl 

/X  Sq  4“ 

1 

u2 

V2  + 

1 

X2 

V2 

1 

yz 

1 

X3 

V3 

(3,9) 
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As 

1  xj  +  ui  yi  +  vi 

1  x2  +  u2  V2  +  v2 

1  X3  +  u3  y3  +  v3 

is  first  order  infinite  small, 


/  = 

Considering  only  the  first  order  small, 


And 


A  =  So  +  ( (m  -  1/3)  (ij -12))  CO 

+  ((»l-»)  (*2  -*l))  • 

(“j)  =P’(*i,!/,)]{D}, 

(“’)  =  [T(*2,j,2)]{Z>}, 

(“’)  =  [T(x3,!o)]{Z)}, 


(3.11) 


(3.12) 
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then 


(3.13) 


A  =  So  +  ( (s/2  -  2/3)  («3  -*2))  [T{xi,yi)]{D}, 

+  ((j/3  -yi)  (xi  —  £3)  )  [T(x2, 3/2 )] {-^} > 

+  ( (2/1  —  2/2)  (X2  -  xi))[T(x3,y3)]{D}. 

Denote 

w=7px*1.».)r(s:s). 

{G}=i(T(x2,!,2)r  (»:») 
+i[T(x3.wr(«:fi), 

and 

w}  =  ipj(«.,».)]r  (»:*). 

{G.}=i[Ti(x2^r(^:^) 

+  )[T,(x3,y,)]T("‘_^J. 

SO 


/*\ 

(°'\ 

#2 

G2 

w  = 

#3 

{G}  = 

G3 

\Hn) 

\6J 

d={ff}T{fl}  +  {G}T{D}  +  ^. 
The  potential  energy  of  the  normal  spring  is 


(3.14) 


(3.15) 


(3.16) 
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n,  =  Id2 


=  t;({H)T{D}  +  {G)T{D)  +  j- 

{D}t{H}{H}t{D}  +  { D}t{G}{G}t{D }  +  2  {D}t{H}{G}t{D} 


2n 


(3.17) 


Thus, 


p{H}{H}T 

p{H}{G}t 

p{G){H}t 

p{G}{G}r 


are  the  normal  spring  matrices  of  [K\ 


~P 

-V 


are  the  load  matrices  of  [F]. 
Then 


p{Hr}{H3}T 

[A  r4] , 

r,  s  =  1, 2, 3, . . .  ,n 

P{Hr}{Gs}T 

[Krs], 

r,  s  =  1, 2, 3, . . . ,  n 

p{Gr}{H3}T 

— y 

[Krs], 

r,  s  =  1, 2, 3, . . . ,  n 

p{Gr}{Gs}T 

[Kr.], 

r,  s  =  1, 2, 3, . . . ,  n 

-p(t)w 

{ Fr }, 

r  =  1, 2, 3, . . . ,  n 

-'’(f)  <G-} 

-> 

{Fr}, 

r  =  1,  2,  3, . . . ,  n 

(3.18) 
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3.5  Shear  Contact  Matrix  for  General  Covers 


In  Figure  3.8,  point  (x0,  yo)  is  on  the  edge  P2P3.  Point  (x0,  y0)  is  also  the  assumed 
contact  point  of  vertex  Pi .  The  shear  spring  on  the  direction  P2P3  connects  vertices  Pi  and 

Po- 


Let 

l  =  v(x2  +  U2-  X3  -  U3)2  +  (2/2  +  t>2  -  2/3  -  t»3)2, 


the  shear  displacement  of  Po  and  Pi  along  line  P2P3  is 
d=jPj^-Pj%  (3.19) 

- )  (v.  +«!)-(».+ «o) )  ( 1  +  $ )  • 


The  step  displacements 

(u{,vi)  i  =  1,2,3 

is  small  as  the  results  of  small  time  step.  The  contact  distance  j  is  small  from  the  definition 
of  the  contacts.  Then  the  projection  of  the  displacements  (u0,  i>o)  and  (ui,v\)  of  points  Po 
and  Pi  on  vector  P2P3 


(x3  +  u3  -  x2  -  u2 
(  X3  +  u3  -  x2  -  u2 


Vz  +  V3  -  y2  -  v2  ) 
yz+v3-y2-v2) 


ai'e  small.  And 

is  second  order  small. 

The  contact  distance  is  small  from  the  definition  of  the  contacts.  Then  P0P1  is  small, 
and 


(  u3  -  u2  v3-v2) 


Uq 

Vo 


{  X3  -  x2  \ 

\y3-y2 ) 


So  =  (  Xi  -  x0  2/1  -  yo  ) 
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is  small. 


( *1  -  x0 


yi-yo) 


f  u3  -  u2 

\Vs  -V2 


is  second  order  small. 


y/{X2  +  U2  -  Z3  -  U3)2  +  (?/2  H-  f 2  —  J/3  —  W3)2 

.  .  .  .  ,  ,  \  \  f  (x3  +  W3)  —  (#2  +  u2) 

((*1  +«i)-(*o  +  uo)  (yi+ui)-(yo+uo))^  (y3  +  V3)_(y2+t?2) 

d  ~ 

_ 1 _ 

\/(x2  -  Z3)2  +  (2/2  -  J/3)2 

( (an  +  tii)  —  (*o  +  wo)  (yi  +  »i)  -  {yo  +  «o) )  ( ^  +  v*) 

Therefore  the  formula  of  l  can  be  simplified  as 


l  =  %/(x2  -  x3)2  +  (j/2  ~  2/3  )2 


(3.20) 


d 


So 

l 

So 

l 


+  T  (  x3  -  x2 
1 , 

+  -  (  x3  -  x2 
1 . 

+  T  (  X2  -  X3 


V3  ~  V2  ) 


f  U\  —  Uo 
\V1-V0 


y*  -  V2 ) 


V2  -  V 3 


(3.21) 


Denote 

w-ipx.„»r  <3-22> 

<3-23> 
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and 


W}=i[T,(x„W)r(;-^), 

{ei}=ip«^)r(2:2). 


(H  i\ 

(G  i\ 

h2 

{H}  = 

h3 

{G}  = 

\Hn) 

\Gj 

Then  (3.21)  becomes 

d=  {H}t{D}  +  {G}t{D}  +  ^.  (3.24) 


The  potential  energy  of  the  shear  spring  is 


n,  =  ?d2 


_  P 
2 

=  t 
2 

+  2 


({H}T{D}  +  {G}t{D}  +  ^ 


{D}t{H}{H}t{D}  +  {D}t{G}{G}t{D}  +  2{D}t[H}{G}t{D} 


s0 


So 


{Dy{H}+2[^){D}I{G}  + 


So 


(3.25) 


Thus, 

p{H}{G}T 

p{G){H}T 

P{G}{G}T 
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are  the  shear  spring  matrices  of  [K] 


are  the  load  matrices  of  [F]  from  shear  spring. 


p{Hr}{H,}T 

[K„l 

p{Hr}{G,}T 

— y 

[Krsl 

p{Gr}{H,}T 

[Kr.], 

— y 

[Kr3], 

— y 

{ Fr }, 

-pfyW) 

-> 

{Fr}, 

r,s  =  1,2,3, ...  ,n 
r,s  =  1,2, 3, . . .  ,n 
r,s  =  1,2,3,. ..  ,n 
r,s  =  1,2,3,. . .  ,n 

r  =  1,2,3,. ..,n 
r  =  1, 2, 3, . . . ,  n  (3.26) 


3.6  Friction  Force  Matrix  for  General  Covers 

Friction  forces  are  treated  as  loading  forces  in  forward  analysis,  therefore  the  coef¬ 
ficient  matrix  of  the  equations  will  still  be  symmetric. 

When  Coulomb’s  law  allows  sliding  between  two  sides  of  boundary  contacts,  there 
exist  friction  forces  in  two  sliding  sides  if  the  friction  angle  4>  is  not  zero. 

As  shown  in  Figure  3.8,  Pi  is  on  one  side  of  a  contact  and  P2,  P3,  Po  tiro  on  the  other 
side  of  the  same  contact.  The  friction  force  is  calculated  from  the  normal  contact  compres¬ 
sive  force  and  the  direction  of  the  friction  force  is  depending  on  the  movement  of  Pi  relative 
to  P0  in  the  direction  from  P2  to  P3.  Let  p  be  the  stifftiess  of  normal  contact  spring,  then 
the  friction  force 

F  =  p  ■  d  ■  s  ■  tan (</>),  (3.27) 


48 


FIGURE  3.8  Position  of  shear  spring 


where 

d  is  the  normal  penetration  distance, 
tan (</>)  is  the  friction  coefficient, 

s  =  sgn  (movement  of  P\  relative  to  P0  in  the  direction  from  P2  to  P3) ,  and 

(  1,  if  x  >  0, 
sgn(x)  =  <  0,  if  x  =  0, 
l  —1,  if  x  <  0. 

The  friction  force  T  is  along  the  direction 


y((x3-x2)  (y3-y2)), 

and 

l  =  v/(2;2  -  z3)2  +  (y2  -  y3)2 
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is  the  length  of  line  P2P3- 


Then  the  potential  energy  of  friction  force  T  at  P\  side  is 


Denote 


tt  f  (  s(x3~xi 

11/  =  —  (  Ui  Vi  ) 

1  l  V  2/3  —  2/2 


(3.28) 


(3.29) 


w}=fps(...»  or 


m= 


Tlien 


■T{H) 


/{Hi} 
{Hi} 
-T  {Hi} 


(3.30) 


{H„} 


is  the  loading  matrix  for  point  P\\  and 
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T{Hr]  ->  {Fr} 


r  =  1,2,3, ...  ,n 


Similarly,  the  potential  energy  of  friction  force  T  at  P0  side  is 


*>(&:$) 


Denote 


{G>  =  j 


and 


SO 


/Gl\ 

<?2 

g3 


Then 


+^-{G}  =  +.F 


{^1}  \ 

{C2} 

{^3} 


'  {Gn  }  / 


(3.31) 


(3.32) 


(3.33) 
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is  the  loading  matrix  for  point  Pq-,  and 


+P{Gr) 


{Fr},  r  = 
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Chapter  4 

Finite  Element  Covers  of  Manifold  Method 


4.1  Finite  Covers  Formed  by  Finite  Element  Nodes 
and  Physical  Boundaries 

The  manifold  method  can  perform  the  computations  of  finite  element  method  for 
continuous  material.  After  transferring  the  finite  element  mesh  to  finite  covers  of  manifold 
method,  the  joint  lines  can  be  input  as  new  physical  mesh,  then  the  same  finite  element  mesh 
can  compute  joints  in  the  same  material  volume. 

The  finite  element  meshes  can  be  used  to  define  finite  covers  for  manifold  method. 
Considering  any  node,  all  elements  having  this  node  form  a  mathematical  cover  (called 
“star”  in  algebraic  topology). 

In  Figure  4.1  and  4.2,  the  mathematical  cover  V5  of  node  5  has  three  elements  2  4  5, 
2  5  3  and  3  5  6.  The  mathematical  cover  V\  of  node  1  has  only  one  element  12  3  which  is 
the  only  element  having  node  1. 

mathematical  covers  of  Figure  4. 1  and  Figure  4.2 


node 

element 

element 

element 

element 

element 

element 

1 

1,2,3 

2 

1,2,3 

2,4,5 

2,5,3 

3 

1,2,3 

2,5,3 

3,5,6 

4 

2,4,5 

5 

2,4,5 

2,5,3 

3,5,6 

6 

3,5,6 

Figure  4.3  shows  the  mathematical  covers  of  nodes  1,2,4, 5,  mathematical  covers 
of  all  nodes  are  listed  in  the  table. 
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mathematical  covers  of  Figure  4.3 


node 

element 

element 

element 

element 

element 

element 

1 

1,2,5 

1,5,4 

2 

2,3,6 

2,6,5 

2,5,1 

3 

2,3,6 

3,7,6 

4 

1,5,4 

4,5,8 

5 

1,5,4 

1,2,3 

2,6,5 

5,6,9 

5,9,8 

4,5,8 

6 

2,6,5 

2,3,6 

3,7,6 

6,7,10 

6,10,9 

6,9,5 

7 

3,7,6 

6,7,10 

8 

4,5,8 

5,9,8 

9 

5,9,8 

5,6,9 

6,10,9 

10 

6,7,10 

6,10,9 

Any  original  finite  element  is  the  common  area  of  the  mathematical  cover  of  its 
nodes.  In  Figure  4.1  and  Figure  4.2,  mathematical  cover  V5  of  node  5  is  the  area  defined 
by  the  polygon  2  4  5  6  3;  mathematical  cover  V2  of  node  2  is  the  area  defined  by  the  poly¬ 
gon  12  4  5  3;  mathematical  cover  V3  of  node  3  is  the  area  defined  by  the  polygon  1  2563. 
The  common  part  of  mathematical  covers  V5,  V2  and  V3  are  original  finite  element  5  2  3. 

The  physical  mesh  of  Figure  4.1  and  4.2  are  the  thick  lines,  which  are  the  bound¬ 
aries  and  the  fractures  of  the  material  volume.  The  physical  covers  or  covers  are  defined  as 
following: 

[  1]  The  region  of  physical  cover  is  the  materials  contained  in  the  mathematical  cover,  or 
mathematically  speaking  the  intersection  of  the  mathematical  cover  and  the  material 
field. 

[2]  If  the  material  boundaries,  block  boundaries  or  fractures  divided  the  mathematical 
cover  to  totally  isolated  regions,  each  region  is  a  physical  cover.  Therefore  the  phys¬ 
ical  covers  are  the  subdivision  of  the  mathematical  covers. 

The  mathematical  covers  and  divided  physical  covers  of  Figure  4.1, 4.2  and  4.3  are 
listed  by  the  tables. 

the  mathematical  covers  and  physical  covers  of  Figure  4.1 
math  cover  physical  cover  physical  cover  physical  cover  physical  cover 
Vi  li 


57 


y2  2X  22 

V3  3i  32 

V4  4i  42 

V5  5i  52 

V6  61  62 


the  mathematical  covers  and  physical  covers  of  Figure  4.2 
cover  physical  cover  physical  cover  physical  cover  physical  cover 
li 

2i  22 

3i 

4i  42 

5i 

61 

the  mathematical  covers  and  physical  covers  of  Figure  4.3 


math  cover 

physical  cover 

physical  cover 

physical  cover  physical  cover 

Vi 

li 

I2 

V2 

2i 

22 

v3 

3i 

v4 

4i 

42 

V5 

5i 

52 

53 

V6 

61 

62 

63 

v6 

7i 

v6 

81 

82 

v6 

9i 

92 

93 

V6 

10i 

4.2  Elements  as  The  Common  Part  of  Node  Covers 

The  elements  of  the  manifold  are  the  common  regions  or  the  intersections  of  the 
physical  covers.  Each  point  inside  the  material  boundary  lies  in  a  “element”  which  is  a  com¬ 
mon  part  of  exactly  three  physical  covers  in  the  triangle  finite  element  case.  The  relations 
of  the  original  finite  element  terminology  and  its  manifold  generalization  are 
dimensions  finite  element  method  manifold  method 

0-d  to  2-d  node  physical  cover 


math 

Vi 

V2 

V3 

V4 

v5 

V6 
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1- dto  1-d 

2- d  to  0-d 


edge 

element 


two  cover  intersection 
intersection  of  covers 


Under  manifold  method,  the  “elements”  and  “nodes”  here  are  the  extensions  of  their 
FEM  counterparts.  Using  the  new  nodes  and  elements,  the  joints  can  open  and  slide,  the 
blocks  can  move  away  and  the  continuous  area  of  the  material  body  can  still  be  connected. 

The  proof  of  these  important  conclusions  come  directly  from  the  definition  of  the 
finite  cover  systems  and  the  cover  functions  and  global  displacement  functions  of  the  man¬ 
ifold  method. 

A  original  finite  element  can  be  divided  to  several  elements  of  manifold  method  by 
the  block  boundaries,  joints  or  fractures.  Before  the  deformation,  the  nodes  share  the  same 
position.  It  can  be  understood  as  many  layers  divided  by  discontinuities  on  the  original  sim¬ 
ple  finite  element  mesh. 

The  manifold  elements  of  Figure  4.1,  Figure  4.2  and  Figure  4.3  are  listed  bellow. 

Physical  covers  or  nodes  of  elements  of  Figure  4. 1 


element  number 

physical  cover 

physical  cover 

physical  cover 

1,2,3 

li,  2i,  3i 

2,4,5 

22, 4i , 5i 

2i,42,52 

2,5,3 

2i ,  52 , 3i 

22,5i,32 

3,5,6 

3i,52,62 

32,5i,6i 

Physical  < 

covers  or  nodes  of  elements  of  Figure  4.2 

element  number 

physical  cover 

physical  cover 

physical  cover 

1,2,3 

li,2i,3i 

2,4,5 

22,4i,5i 

2j,42,5i 

2,5,3 

2i,51,31 

22, 5i , 3i 

3,5,6 

3i,5i,6j 

Physical  covers  or  nodes  of  elements  of  Figure  4.3 

original  element 

physical  cover 

physical  cover 

physical  cover 

1,5,4 

11,51,4! 

1,2,5 

li,  2i ,  5i 

12,22,52 

2,6,5 

2i,6i,5j 

22, 62, 52 

2i ,  63 , 5i 

2,3,6 

2i,  3i,6i 
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3,7,6 

3i,  ?i ,  6i 

4,5,8 

4i,5i,82 

42, 53 , 8i 

5,9,8 

5i  ,  92, 82 

53, 9i ,  8i 

5,6,9 

53?  6i , 9i 

5i,63,92 

5,6,9 

51,61,94 

6,10,9 

6i,10i,94 

6i,10i,9i 

6,7,10 

6i,7i,10i 

In  Figure  4.1,  the  only  joint  inside  the  material  divided  completely  the  material  to 
two  disconnected  parts,  any  manifold  element  divided  by  this  joint  have  completely  differ¬ 
ent  nodes  or  physical  cover  numbers.  Therefore  the  two  manifold  elements  are  free  to  move 
independently. 


element  separation  along  joint  in  Figure  4.1 
original  element  manifold  element 

2.4.5  22,4i,5i 

2,5,3  2i,52,3i 

3.5.6  3i,52,62 


manifold  element 

2i ,  42, 52 

22. 51 ,  32 

32. 51 ,  6i 


In  Figure  4.2,  the  only  joint  inside  the  material  divided  partly  the  material,  any  man¬ 
ifold  element  divided  by  this  joint  have  partly  different  nodes  or  physical  cover  numbers. 
Therefore  the  two  manifold  elements  can  have  different  movements. 


element  separation  along  joint  in  Figure  4.2 
original  element  manifold  element  manifold  element 

2.4.5  22,4i,5i  2i,42,5i 

2,5,3  21,5i,31  22,5i,3i 

3.5.6  3i,5i,6i 

If  two  manifold  elements  share  a  edge  of  original  finite  element  mesh,  as  this  edge 
is  not  a  joint,  the  two  manifold  elements  have  the  same  nodes  on  the  edge: 


element  connection  along  edges  in  Figure  4.1 


nodes  of  edge 

manifold  element 

manifold  element 

2,3 

li,2i,3i 

2i,52,3i 

2,5 

2i,42,52 

2i,  52, 3i 

2,5 

22,4i,5i 

22,5i,32 
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3,5 

2i,  52, 3i 

3i,  52, 62 

3,5 

22 , 5i ,  32 

32,5i,6i 

element  connection  along  edges  in  Figure  4.2 

nodes 

of  edge  manifold  element 

manifold  element 

2,3 

1  i  )  2i ,  3i 

2i,5i,3i 

2,5 

2i,42,5i 

CO 

iO 

<N 

2,5 

22,4i,5i 

22, 5i ,  3i 

3,5 

2i,5i,3i 

3i,5i,6i 

3,5 

22, 5i,  3i 

3i,  5j , 61 

Using  the  manifold  definition  of  nodes  and  elements,  the  following  important  con¬ 
clusions  can  be  proved  and  also  can  be  seen  directly  from  Figure  4.1,  4.2  and  4.3: 

[  1]  the  elements  are  irregularly  shaped; 

[2]  each  element  has  three  physical  cover  numbers; 

[3]  these  three  covers  have  one  element  as  their  common  area; 

[4]  the  three  covers  can  be  seen  as  three  “nodes”  of  the  element; 

[5]  the  adjacent  elements  have  the  same  nodes  along  the  common  edge; 

[6]  two  elements  divided  by  fractures  or  boundaries  have  different  nodes. 


4.3  Cover  Functions  and  Weight  Functions  of  Finite  Element  Mesh 

The  displacement  functions  are  independent  from  the  material  boundary.  If  the  ma¬ 
terial  occupies  only  part  of  the  element,  the  displacement  functions  are  still  the  same.  For  a 
triangular  element,  there  are  three  covers  containing  this  element  corresponding  three  nodes. 
Therefore  each  element  is  the  common  region  of  three  covers  of  its  three  nodes. 

For  a  element,  the  weight  functions  are  computed,  denote  ?,•  :  ( )  the  coordi¬ 
nates  of  nodes  i  =  1, 2, 3,  and  the  related  nodal  displacements  as  follows: 
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coordinates 

i\  ■  (a?i,yi) 

*2  :  (*c2 1  Vi ) 
*3  :  (®3,y3) 


displacements 

(u2,^2), 

(«3,«3)- 


The  displacement  field  can  be  described  as: 


u  =  ai  T  &ix  +  ciy, 
t;  =  a2  +  b2x  +  C2  y , 


(4.1) 


The  nodal  displacements  are 


«i  =  <*i  +  hxi  +  cij/i, 
U2  =  «i  +  a:  2  +  cij/2, 

^3  =  ^1+  b  1X3  +  C1J/3, 
Ul  =  02  +  ^2^1  +  C2  J/l , 
U2  —  02  +  b2X2  +  C2y2, 
V3  =  02  +  biXz  +  C2y  3, 

and  by  matrix  notation  equation  (4.2)  becomes 


and 


/ 1  arj  yi  \  /ai  \ 

I  1  a:2  J/2  J  I  I  » 

\  1  ar3  y3  /  \  ci  / 


Then 


/  t>i  \  /  1  a:i  yi  \  / a2  \ 

I  u2  1  =  I  1  a;2  j/2  I  I  &2  1 

\t>3  /  \  1  x3  y3  /  \  c2  / 


/  ai  \  /  1  Xi  yi  \  /ui  \ 

I  61  =  I  1  x2  y2  j  I  u2  )  , 

\cx  /  \  1  x3  y3  /  \U3  / 


(4.2) 


(4.3) 


(4.4) 
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5 


and 


thus. 


and 


Denote 


and 


1  xi  yi 
u-  (l  x  y)  |  1  x2  y2 
1  xz  yz 

1  *i  2/i 
v  =  (1  x  y )  (  1  x2  y2 
1  xz  yz 


1  x1  yx 

( h  h  fz )  =  ( 1  x  y )  [  1  x2  y2 

1  xz  yz 


fn  fn  fiz\  /l  xx  yx 

/21  /22  /23  I  =  I  1  x2  y2 

Jzi  fz2  fzz  /  \1  X3  2/3 


and  the  determinant 


A  = 


1  x2  y j 

1  x2  2/2 

1  x3  yz 


-1 


(4.5) 


(4.6) 


is  two  times  the  area  of  the  element. 


then 


fll  fl2  flZ 

/21  /22  /23 

/31  fz2  fzz 


( 


-1 


1 

A 


+ 


+ 


x2 

V2 

_ 

X\ 

yi 

+ 

X\ 

2/1 

\ 

xz 

V3 

Xz 

2/3 

X2 

V2 

1 

1 

2/2 

2/3 

+ 

1 

1 

yi 

V3 

— 

1 

1 

2/i 

2/2 

1 

1 

x2 

xz 

— 

1 

1 

CO  (-1 

+ 

1 

1 

X! 

^2 

/ 

l  ( x2yz  x3y2  x3yi  -  xiyz  xxy2  -  x2t/i 

^  I  V2~yz  yz  ~  2/1  2/i  -  2/2 

\  Xz  —  X2  X’!  —  Xz  X2  —  XX 


(4.7) 
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(  fll  fl2  fl3 

(/l  h  /*)  =  (!  *  y)(/21  /22  /23 

V/31  fz2  fs3 


and 


//i\  /  /11  +  /21*  +  fsiy  \ 

j  /2  I  —  J  /l2  +  /22^  +  /32?/  I 

V/3/  \/l3  +  /23^  +  /33j/  / 


«  =  (/l 


/2 


h ) 


v  =  {h  h 


h  ) 


(4.8) 


(4.9) 


In  the  finite  element,  the  cover  functions  are  the  constant  function  u,,  Vi  over  the 
whole  cover.  The  cover  functions  (u,(x,  y),Vi(x,  y ))  defined  on  physical  cover  Ui 


f  Ui(x,y) 
1  Vi(x,y) 


(x,y)  G 


(4.10) 


The  weight  function  w, ;(x,  y)  is  the  shape  function  of  finite  element  method. 


Wi(x,y)  =  fi(x,y)  i  =  1,2, ...  ,n 

The  meaning  of  the  weight  functions  tc,  (x,  y)  is  weighted  average,  which  is  to  take  a 
percentage  from  each  node  displacements  u{,  v, .  For  each  element,  the  summation  of  three 
weight  function  of  three  nodes  are  1. 


»,-(*,!()  =  1.  (4.11) 

(x,y)€Uj 

Equation  (4.11)  can  be  derived  from  equation  (4.7), 
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fn  +  /12  +  /13  =  1 
/21  +  y*22  +  /23  =  0 
/31  +  fz2  +  /33  =  0 


wi(x,y )  +  w2(x,y)  +  w3(x,y) 
=fi(x,y )  +  f2(x,y)  +  f3(x,y)  =  1 


4.4  Finite  Element  Global  Functions  for 
Continuous  and  Discontinuous  Materials 

Using  the  weight  functions  to;( x,  y)  a  global  function  on  the  whole  physical  cover 
system  is  defined  from  the  cover  functions 


f  u(x,y) )  =  (T:=1m(x,y)Ui\ 

\u(x,x/)J  \E"=i  Wi(x,y)vi  J 


(4.12) 


/  u(x,  y)  \  _  (  E”=  1  fi(x » 2/)wt  \ 
\u(.T,y)j  vEr^i/ii^yK  / 


(4.13) 


Then 

fUl\ 

vi 

f  u(x,  y)  1  /  f\  0  f2  0  /3  0  \  u2 

\  v(x,  y)  J  \  o  f!  o  f2  o  f3  J  v2 

U3 

\V3  J 

=  [Te]Pe},  (4.14) 


where 


and 

[Ti] 


{Di}\ 

[Te]  =  ([Tx]  [T2]  [T3]),  {De}=  \{D2}  , 

VPs}/ 


( fi(x,y)  0  \ 

V  0  fi(x,y)  J  ’ 


i  =  1,2,3. 


/  h{x,y)\  (  fn  +  hix  +  /3iy  \ 

f2(x,y )  =  fn  +  fi2X  +  f32y 

\f3(x,y)J  \/i3  T  f23x  -T  fzzy  ) 


65 


The  cover  functions  of  the  triangular  finite  element  method  are  constants,  the  global 
functions  are  linear  functions  of  the  coordinates  (x,  y). 

For  manifold  method,  the  cover  function  can  be  normal  two  dimensional  series, 
where  the  unknowns  are  the  coefficients  of  the  series. 


f  Ui(x,y)l  =  /£7=i«^;(*,y)\ 

1  «»•(*,  y)  j  V  £7=  i  vas](x>  y) ) 


(4.15) 


The  global  displacement  function  is  the  combination  of  locally  defined  series. 


f  u(x,y) 1 
1  v(x,y)  f 

(  ELi  wi(x^y)u*(x^y)\ 

/£"=i  £7=1  utjsj(x,y)wi(x,y)\ 
V  £"=i  £y£  vljSj(x,y)wi(x,y)  ) 


The  submatrix  form  of  the  global  displacement  function  is 


f  u{x,y) 
\  v(x,y) 


=  [Te}{De}, 


(4.17) 


where 


{Di}\ 

[T.]  =  ([T1]  [T2]  [T3] ) ,  {D.}=  p2}  , 

\{Dz}J 


and 


mi  = 


66 


(fa  0 

f,  0  /„ 


/,  2  0 

o  1,1 


fiZ  0  .  fnm 

0  fiZ  .  0 


(  ^ 
U,1 
Ui  2 

Vi2 

Ui3 

ViZ 


i  =  1,2,3. 


I  v, im  I 

\Vim  ) 


fij{xi  y)  =  * j(x,y)wt(x,y)  i  =  1,2,3, 


j  =  1,2,3,...,  m. 


4.5  Finite  Element  Equilibrium  Equations  Based  on  Manifold  Method 

The  total  potential  energy  II  is  the  summation  over  all  the  potential  energy  sources: 
individual  stresses  and  forces.  In  the  following,  the  potential  energy  of  each  force  or  stress 
and  their  differentiations  are  computed  separately. 

[1]  the  strain  potential  energy  II £  produces  the  stiffness  matrix, 

[2]  the  potential  energy  Ua  of  initial  stresses  produces  the  initial  stress  matrix, 

[3]  the  potential  energy  IIP  of  point  load  produces  the  point  load  matrix, 

[4]  the  potential  energy  !!„,  of  body  load  produces  the  body  load  matrix, 

[5]  the  potential  energy  11,  of  inertia  produces  mass  matrix, 

[6]  the  strain  potential  energy  II3  of  contact  springs  produces  contact  matrices, 

[7]  the  potential  energy  11/  of  friction  forces. 

There  are  n  physical  covers  or  n  nodes.  For  the  two  dimensional  triangle  elements 
j,  there  are  three  physical  covers  or  nodes  per  element  ji,j2,jz-  Each  physical  cover  or 
node  i  has  two  unknowns  (tt^,  v,) 
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Pi]  = 


dn 

di2 


Therefore  the  unknown  [D,]  of  cover  i  is  a  2  x  1  submatrix.  The  coefficient  matrix 
should  formed  by  the  2x2  submatrices  [Kij]  in  the  equation  (4.18).  The  force  term  Ft  on 
cover  i  is  also  a  2  x  1  submatrix  due  to  the  dimension  of  Dz 

Adding  all  potential  energy  together,  the  total  potential  energy  has  the  form 


(Kn 

Ku 

Ku 

Kin  \ 

(Dl\ 

n  =  l(Bf  Dj  Dj  . 

f  I<21 

K2  2 

k23 

...  I<2n 

d2 

■  DTn) 

Kzi 

k32 

K33 

...  Kin 

D3 

\Knl 

Kn2 

Kn3 

...  Knn) 

\DnJ 

+  {DJ  Dj  Dj  ...  D 


'S' 

F3 

\fuJ 


c 


(4.18) 


From  the  formulation  of  II,  the  formula  (1.12)  can  be  written  as  a  symmetric  repre¬ 
sentation, 

Kij  =  Kji 


For  cover  i,  equations 


an 

dvi 


(4.19) 


represent  the  equilibrium  of  all  the  loads  and  contact  forces  acting  on  cover  i  or  node  i  along 
x  and  y  directions  respectively. 


The  differentiations 


a2n 


r,s  =  1,2, 


ddirdd __ 

are  the  coefficients  of  unknowns  djs  of  the  equilibrium  equation  (4.18)  for  variable  d 


(4.20) 


The  differentiations  (4.20)  of  the  total  potential  energy  n,  produces  n  equations  of 
submatrices  as  there  are  n  physical  covers  or  nodes  in  the  manifold.  The  simultaneous  equi¬ 
librium  equations  have  the  same  coefficient  matrix  as  the  quadratic  form  of  II  in  (4.18): 
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Because  each  node  or  cover  has  two  degrees  of  freedom  in  2-d  FEM  manifold,  each  element 
[Kij]  in  the  coefficient  matrix  given  by  equation  (4.21)  is  a  2  x  2  submatrix.  {£>*}  and  {F,-} 
are  the  2  x  1  submatrices. 


In  case  the  cover  functions  are  series  with  m  unknown  coefficients,  the  displacement 
function  ( Ui(x,y),Vi(x,y ))  have  2 m  unknowns.  Therefore  each  submatrices  Kij  in  the 
coefficient  matrix  given  by  equation  (4.21)  is  then  a  2m  x  2m  submatrix.  Dt  and  F,  are 
2m  x  1  submatrices,  where  Dt  represents  the  displacement  variables  of  physical  cover  i. 


69 


70 


Chapter  5 

Element  Matrices  of  Finite  Element  Covers 


5.1  Stiffness  Matrix  for  Finite  Element  Covers 

For  FEM,  the  integration  domains  of  the  stiffness  matrices  are  whole  elements  with 
standard  boundaries.  For  manifold  method,  the  integration  domains  of  the  stiffness  matrices 
are  the  manifold  elements,  which  can  be  part  of  the  elements. 

Same  as  FEM  method,  the  relationship  between  stress  and  strain,  is  given  by 


E 

~  1  -  v2 


V 

1 

0 


/  €x  \ 

~  [E]  (  ey  I  5 
\lxy  / 


(5.1) 


where 

A  ,  °  \ 

I  V  1  o  , 

\0  0  ^ J 


and  E,  v  are  the  Young’s  modulus  and  Poisson’s  ratio  respectively.  And 


(“)=PV]{D.} 

=  ([r,)  [t2]  [rs])  {d,}  , 

\m) 


where 


[T,  1  = 


{B.}  = 


1 


Ui 

Vi 
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and 


/  fi  \  //n  +  fi\x  +  hiy\ 

I  /2  I  =  I  fl2  +  fl2x  +  fziy  1 
V/3/  \  /l3  +  f23x  +  fzzy  ) 


Then 


0 

dx 

0 

a/» 

dx 

0 

0 

0 

dh 

0 

“ “ 

dy 

dy 

dy 

dx 

hi 

0 

f22 

0 

/23 

0  \ 

=  0 

hi 

0 

/32 

0 

/33 

V/31 

hi 

/32 

/22 

/3  3 

/23  / 

/tij\ 

)V1 
U2 
V2 

U3 

\  V3  / 
/til\ 

Vi 

u2 

v2 

u3 

\v3  J 


Let 

[B,]  =  ([B1]  [52]  [£3]), 


where 

//«  °\ 

[Bi]=  0  /3t  I  ,  *  =  1,2,3. 

\/3i  fli  ) 


Then 


=  [Be]{De}  =  ([B1] 


m 


{Di}\ 

[£3])  {£2} 

\{D3}J 


(5.2) 


(5.3) 


(5.4) 


(5.5) 


The  strain  energy  ne  done  by  the  elastic  stresses  of  element  e  is 


ne  =  JJ  ~(ex <jx  +  eyay  +  7*3, TXJ,)  dx  dy , 


where  the  integration  is  over  the  entire  material  area  A  in  that  element.  Then 


(5.6) 
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ne  —  J J*  2  (  ^ y  ^Yxy  )  |  J  ^x  dy 

^  \  TXy  J 

=  \  jJ^{De}T[Be]T[E}[Be}{De}dxdy 
=  \{De}T  JJjBe}T{E][Be]dxdy  {De} 
=  \{De}T(se[Be)T[E][Be\){De}, 


where  Se  is  the  area  of  that  element. 


Therefore, 

f[Bi}T 

Se[Be}T[E}[Be}  =  Se  [. B2]t 

\mT 

is  the  element  stiffness  matrix. 

Then 

Se[Br]T[E}[B3}  -+  [Kl{r)i(s)},  r,s  =  1,2,3, 


(5.7) 


(5.8) 


and 


m 


i  =  1,2,3, 


where 


( *1, 

e  =  i. 

first  node  of  the  element 

i(£)  =  s  *2, 

i  =  % 

second  node  of  the  element 

l  *3, 

£  =  3, 

third  node  of  the  element 

5.2  Initial  Stress  Matrix  for  Finite  Element  Covers 

Following  the  time  sequence,  the  manifold  method  computes  step  by  step.  The  com¬ 
puted  stresses  of  previous  time  step  will  be  the  initial  stresses  of  the  next  time  step.  There¬ 
fore  the  initial  stresses  are  essential  for  manifold  computation. 
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For  the  element  e,  the  potential  energy  of  the  initial  constant  stresses 


{*«}  =  (  *2  <  r°xy) 


is 


n„  =  Jj ' 


+  lxyTxy)dx  dy 


-ywy 


f  ax\ 

ey  Ixy )  e°y  dxdy 

\  d  I 

\  xy  / 


~  JJ  {De}T[Be]T{cr°}dxdy 

=  S*{De}T[Be]T{a°e}, 


where  Se  is  the  area  of  that  element.  Therefore, 


(5.9) 


is  the  load  matrix. 


/[^]T\  /*s\ 

mT 

\[*«]T/  \r%) 


Then 


-Se[Br]T 


{ ( r )  }  t  r  —  1,2,3. 


(5.10) 


where 


f  *1, 

t  =  1, 

first  node  of  the  element 

i(t)  =  <  *2, 

1  =  2, 

second  node  of  the  element 

l  *3  > 

e  =  3. 

third  node  of  the  element 

5.3  Point  Loading  Matrix  for  Finite  Element  Covers 

Different  from  ordinary  FEM  method,  a  load  point  can  be  any  point  in  its  element. 
The  point  loading  force 
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acts  on  point  (x,  y)  of  element  e.  And 


(F,  Fy  )T 


( u\  _  ( 

\V  )  \v(x,y)  J 


The  potential  energy  due  to  the  point  loading  is 


lip  =  ~(Fxu  +  Fyv ) 

=  -(u  .)(£) 

=  -{£>«}r[T«(x,!/)]T(^. 


Therefore, 


is  the  load  matrix. 

Then 

PV]r(£')  (F«)’  ’•  =  1.2.3. 

where 


f  *1, 

i  = 

1, 

first  node  of  the  element 

*(^)  =  {  *2, 

£  = 

2, 

second  node  of  the  element 

l  *3, 

1  = 

3, 

third  node  of  the  element 

(5.11) 


(5.12) 


5.4  Body  Loading  Matrix  for  Finite  Element  Covers 

Assuming  that 


(/*  f,f 


75 


is  the  constant  body  force  acting  on  the  material  area  of  element  e. 


The  potential  energy  due  to  the  body  loading  is 


nw  = 


(fxU  +  fyv)  dx  dy 


j  (u  u )  I  *  )  dx  dy 
A  \Jy ) 


-{DeY 


I  [Tey  dx  dy  (  x 
A  J  \Jy 


(5.13) 


Therefore, 


f  [Tey  dxdy  / 
A  J  \Jy 


is  the  body  loading  matrix.  And 


f  ([TjM  dx  J  (f 
A  V  IT.  IT  /  V  h> 


(5.14) 


[T{]t  dx  dy  = 


Si  0 


dx  dy 


'A\0  fi  J  y 

[  (  fu  +  hix  +  fay 


fu  +  hix  +  f3ly 


dx  dy 


fuSe+hiSex  +  f3iSt 


fliSe  +  /2  i$x  +  hiSy 


Then 


JL™ 


dx  dy  = 


St  0  \ 

o  Si  r 


where 


5,-  =  hiSe  +  hiSl  +  hiSh 


se  = 


ce  _ 
— 


ce  _ 


dx  dy , 


x  dx  dy , 


y  dx  dy. 


(5.15) 

(5.16) 


(5.17) 
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Then 


(Sr  0  \(f,\ 

\  0  Sr  )  \fy  ) 


{  F'i  ( r)  }  >  *  —  1,2,3. 


where 


•'(*) 


i'i ,  £  —  1,  first  node  of  the  element 

i2 ,  £  =  2,  second  node  of  the  element 

,  £  =  Z,  third  node  of  the  element 


5.5  Inertia  Force  Matrix  for  Finite  Element  Covers 

The  inertia  force  matrix  is  equivalent  to  the  mass  matrix  of  FEM.  This  matrix  is  the 
most  important  matrix  of  manifold  method.  Giving  a  small  time  steps,  the  inertia  force  ma¬ 
trix  will  control  the  movements  of  all  points  of  the  whole  material  volume.  Choosing  small 
time  steps,  the  discontinuous  contact  computation  will  be  stable. 

Considering  the  current  time  step,  denote 

(u{t)  v(t))T 


as  the  time  dependent  displacements  of  any  point  (x,  y)  of  element  e  and  M  as  the 
mass  per  unit  area.  The  force  of  inertia  per  unit  area  is 


fx 

fy 


,.d 2  (n{t)\  ,.,T,d2{De(t)} 

=  dfl . 


The  potential  energy  of  the  inertia  force  of  element  e  is 


(5.18) 


ni =-//.<"  •>(£)* 


dy 


M  (  u  v)[Tef2{Dd^dxdy. 


(5.19) 


Assume  {£>e(0)}  =  {0}  as  the  element  displacements  at  the  beginning  of  the  time 
step,  {De(A)}  =  {De}  as  the  displacements  at  the  end  of  the  time  step,  A  as  time  interval 
of  this  time  step.  Then 
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P«)  =  {D,(A)}  =  {-0,(0)}  +  Aa{^f(0)}  +  ^a2{g;2(°)} 

3{Pt(0)}  A^(Z?,(0)} 

dt  2  dt2  ’ 


(5.20) 


a2wo)}  _  a(Z)  i  _  imm  =  ±w]_  i 

df2  "A2^  A  dt  A2i  e)  A1  M 


(5.21) 


where 


{K(0)}  = 


a{o,(p)} 

dt  ’ 


(5.22) 


is  the  velocity  at  the  beginning  of  the  time  step.  Then  the  potential  energy  becomes 


n t  =  M{De}T  JJjTe}T[Te}dxdy  |{V«(0)} 


(5.23) 


Therefore, 


-  7/;ri 


[Te]  <ix  dy 


^W}-T<y*(o>} 


is  the  load  matrix.  Thus, 


(5.24) 


dx  dy 


[[  fkrVmi  py  [T3])dxdy 

JJa  \[T*]t) 


[Ke],  (5.25) 


JJjTe}T[Te]dxdy  { Ve(0)} 


f[  [T2V  ([Ti]  [T2]  [T3] )  dx  dy  (K(0)} 

y'/-4  V[r3]r/ 


{Fe}.  (5.26) 
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Then 
2  M 


A2 


JJ[Tr\T[Ta]dxdy 


[Ki  (r)i(s)]>  r,s=  1,2,3.  (5.27) 


2  M 
A 


JJ[Tr]T[T3)dxdy 


where 


WO)}  ->  {Fj(r)},  (rs  (5.28) 


{ v  fo)>  =  —  f 

1  dt \VS(0) 


The  matrix  integration  can  be  computed  as 

[Tr]T[Ta}dxdy 


II, 
-//.( 


US,  0 

'.4  V  0  frf. 


dx  dy 


trs  o 

0  tr3 


where 


=  JJ  frfs  dx  dy 

-  JJ  (fir  +  f2rX  +  f3ry)(fls  +  f2sX  +  hsy)dxdy, 


and 


where 


■flrflsSe 

+  (flrf2s  +  fl»f2r)Sl 
+  (fir  f 3a  +  flsf3r)Sey 

+  f2rf2sSexx 

+  f3rf3sSyy 
+  (/2r/3s  +  f2sf3r)Sexy 

(5.29) 

Ky  =  JJ^xydxdy, 

(5.30) 

>ii=  ff  x2  dx  dy, 

(5.31) 

J  Ja 

'ly  =  JJAy2  dxdy- 

(5.32) 
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As  this  element  of  the  manifold  method  is  a  generally  shaped  polygon, 


Pi  P2  ...  Pm- 1  Pm  Pi,  Pi  =  Pm+1,  P%  =  j,  2/i ) 

are  its  ordered  vertices  rotating  from  axis  x  to  axis  y.  Denote  Pq  =  (x0,  yo)  as  the  arbitrary 
chosen  point,  The  analytical  solutions  of  these  integrations  are  the  following: 


s'=  E 


1  =  1 
<0  m 


1  yo 

1  Xi  yi 

1  Xi+ 1  Vi+i 


Sx  =  —  X^(xo  +  Xi  +  xi+i ), 


Z  —  1 


=  “o"  ^  yi  + 


3  ^ 

5e 


1=1 

m 


xx  g  /  4~  4""  ^z+l  4~  xix0  4”  ^z+1^0  4“  XiZi+i), 

z  =  1 
£e  m 

=  Tj-  + y^  +  ^2+i  +  y*y°  +  Vi+iyo  +  y»y*+i), 


6  2 
S 


t=i 

m 


^xy  ^  ^  "(2xoyo  P  T  2xt_f  iy,+i 

*=l 

P  Xiy0  p  Xi-^-iyo  4*  x0yt  -f-  Xoyi-i-i  T  .Tiyi-f-i  P  ‘^'i-riy*)* 

Then  the  final  formula  is 


(5.33) 

(5.34) 

(5.35) 

(5.36) 

(5.37) 


(5.38) 


2 M  ftra  0 

A2  I  0  Us 


[■^i(r)i(j),,  r,s  —  1,2,3,  (5.39) 


2^Co 


r  =1,2,3, 
s  =  1,2,3. 


(5.40) 


where 


r‘i ,  i  =  1,  first  node  of  the  element 

i(f)  =  •{  i’2,  £  =  2,  second  node  of  the  element 

*3 ,  ^  =  3,  third  node  of  the  element 
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5.6  Fixed  Point  Matrix  for  Finite  Element  Covers 


As  a  boundary  condition,  some  of  the  elements  are  fixed  at  specific  points.  The  con¬ 
straint  can  be  applied  by  using  two  very  stiff  springs.  Assume  the  fixed  point  is  (x,  y)  at 
element  e  with  the  displacements 


/ u(x,y)  \  _  /0\ 
^u(x,y)y  \0  )  ' 


The  computed  displacements  are 

/  \T 
(u  V ) 


There  are  two  springs  which  are  along  the  x  and  y  directions  respectively.  The  stiff¬ 
ness  of  the  springs  is  p.  The  spring  forces  are 


u 

fy 


—pu 

—pv 


The  strain  energy  of  the  spring  is  II/,  then 


n/  =  |  (n2  +  v2) 

=  l(u  v)(U 
2  v  ;  \v 

=  \{Dc)T[Tc\Tm{Dc}. 


(5.41) 


Then 


p[Tey[Te}=p\  [T2 


IT 


[TslT 


pi]  [Ta]  [Ts]) 


is  the  [Ke]  matrix.  Therefore, 


(5.42) 


Pfrfs 


1  0 
0  1 


-» 


[■^i(r)t(s)] 


r,s  =  1,2,3. 
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where 


1  =  1, 

first  node  of  the  element 

*2, 

1  =  2, 

second  node  of  the  element 

*3, 

£  =  3, 

third  node  of  the  element 

Chapter  6 

Entrance  Theory  on  Contacts  of 
Finite  Element  Covers 


6.1  Contacts  and  kinematics  for  Finite  Element  Covers 


As  the  physical  covers  are  formed  by  the  triangle  element  meshes  and  physical 
boundaries,  the  contact  surfaces  are  edges  from 


[1]  block  boundaries, 

[2]  one  of  two  sides  of  joints  in  the  material  field. 


The  vertices  of  the  edges  can  be 

[1]  the  vertices  of  physical  boundaries  or  joints, 

[2]  the  intersection  points  of  triangular  with  block  boundaries  or  joints. 


The  simple  penetration  inequality  in  Chapter  3 


1  xi+ui  y\+v\ 
1  x2  +  u2  J/2  +  v2 
1  a '3  +u3  y3  +  v3 


<  0. 


(6.1) 


is  fundamental  for  contacts. 

[1]  the  inequality  (6.1)  always  correct  when  three  points  Pi ,  P2,  P3  moves  simultane¬ 
ously, 

[2]  the  inequality  (6.1)  is  still  accurate  when  rotation  takes  place, 

[3]  the  assumption  of  small  displacements  (ui,v{)  is  not  needed, 

[4]  no  second  order  infinite  small  is  omitted. 

It  seems,  the  judgment  of  three  point  penetration  is  a  simple  matter.  One  can  find 
several  formulae  for  the  same  judgment.  From  mathematics  the  inequality  (6.1)  is  the  only 


83 


complete  accurate  formula.  Neglecting  second  order  small  of  other  formulae  are 

also  correct  and  accurate  if  the  movements  are  small  enough. 

The  formula  of  contact  distance  (3.3)  of  Chapter  3  is  also  an  accurate  formula. 


1 

1 

1 


*1  +«i  2/i  +  vi 

X2  +  U2  t/2+  V2 
X3  +  ^3  2/3  +  ^3 


(6.2) 


l  =  y/{x2  +  U~2  ~  ~  «3)2  +  (2/2  +  V2  -  2/3  “  U>)2 


However,  in  (6.2),  formula  of  /  is  complex  and  should  be  simplified.  Under  assump¬ 
tion  of  small  pre-existing  distance 


d  = 


1 

y/(x2  -  XzY  +  (2/2  ■  2/3 )2 


1 

Xi 

yi 

1 

X2 

V2 

1 

xz 

Vi 

Distance  formula  can  be  simplified  as 


d  = 


_ 1 _ 

V/(x2  -  x3)2  +  (y2  ^  2/3  )2 


1 

Xi  +  Ui 

2/1  +  Vi 

1 

X2  +  U2 

2/2  +  V2 

1 

X3  +  U3 

2/3  +  ^3 

(6.3) 


In  formula  (6.3),  the  denominator  does  not  include  unknown  variables  (u,,u;). 
Therefore,  formula  (6.3)  is  a  simple  function  with  respect  to  unknown  variable  vt). 


The  inequalities  (6.1)  will  be  a  equation  having  “=”  sign  instead  of  “<”  sign  if  its 
contact  becomes  a  close  contact.  A  stiff  spring  will  be  added  on  the  contact  to  fulfill  the 
equation  (6.1). 


The  friction  law  or  Coloumb’s  law  is  also  inequalities. 


T  <  Ar  tan(0)  +  C 


(6.4) 


The  solution  of  this  inequality  is 
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(  apply  a  shear  stiff  spring  if  IF  <  Af  t&n((j>)  +  C 

|  apply  a  pair  of  friction  forces  ifT  —  A/"tan(<^)  +  C 

In  this  way,  the  inequality  of  friction  law  is  also  transferred  to  equations. 

The  kinematics  of  the  triangle  elements  is  the  same  as  the  general  kinematics  in 
Chapter  3: 

[1]  distance  criteria  of  contacts, 

[2]  angle  criteria  of  contacts, 

[3]  definition  of  entrance  lines, 

[4]  penetration  judgments, 

[5]  penetration  prevention, 

[6]  entrance  position. 


6.2  Normal  Contact  Matrix  for  Finite  Element  Covers 

Assume  Pi  is  a  vertex,  P2P3  is  the  entrance  line  and  ( xk,yk )  and  (ujt,  Vk)  are  the 
coordinates  and  displacement  of  P* ,  k  —  1,2,3  respectively.  If  points  Pi ,  P2  and  P3  rotate 
in  the  same  sense  as  the  rotation  of  ox  to  oy  (see  Figure  3.3  and  3.4),  then  the  distance  d 
from  Pi  to  entrance  line  P2P3  is 


d  — 


A 

7 


1  x\  +  u\  2/1  +  Vi 

1  x2+u2  y2  +  v2  , 

1  x3  +  u3  2/3  +  ^3 


l  =  \/(z2  -  z3)2  +  (2/2  -  y3)2- 


(6.5) 


d  should  be  zero  if  Pi  passed  edge  P2  P3 


l  «i  yi  l  xi  vi 

Ar;S,o+  1  u2  y2  +  1  x2  v2  , 

1  tx3  y3  1  a;  3  u3 


Then 


A  =  So  +  ( (y2  -  y3)  (2:3  -*2)) 

+  ((y3-yi)  (®i-X3)) 


+  ( (yi  —  y2 )  (x2  —  ari) ) 


=  pi(*i,yi  )]{A' }, 


=  [T>(x2,y2)]{T>i}, 
=  Pj  (^3  *  y3 )]  {Dj}, 


A  =  S0  +  ( (y2  -  y3)  (x3  -  x2) )  [Tj(zj ,  yi  )]{A}, 

+  ((ys  -  yi)  (xj  -x3))[Tj(x2,y2 )]{!>>}, 
+  ((yi  -V2)  (x2  -x1))[T7(x3,y3)]{T>J}. 


Denote 


W  =  y[Ti(x1,y1)f 

{G}  =  j[Ti(x2,y2)]: 

+  j[Tj(x3,y3)]2 


T  I  V2-V3 

\X3  -  x2 

it  /  !/3  -  yi 
\Xi  -  x3 

t  (  yi  -  V2 

\x2  -  Xi 


(6.10) 


{Hl}\ 

=  {Hi} 

\{Hs}J 

1  /[^i]T(xi,yi) 

=  7  PiFfo.yi) 

\[^3]T(xi,yi) 


J/2  ~  V3 
X3  -  X2 


86 


\  {^3  }  / 

,  /[Ti]r(x2>!,2) 

=  7  Pi]r(*2,ltt) 

V[r,]T(x2,!,2) 

,  fmT(x„y3) 
+  7  Pi]T(*»,w) 
'  \Pyr(*.,v>) 


T 


2/1  2/2 
£2  “  Zl 


where 


/  ei  \  ^  (  fr{$  1?  2/1 ) 

{Fr}  =  7  V  o 

r^Y  1  _  1  /  fr(x 2,  2/2) 

/  v  0 

1  /  /r(z3,  ^3) 

+  i  v  0 


0  \  (  m  -  V3  \ 

fr(x  1,  yi)y  \x3-x2  J  ' 

0  W  y3  -  y\  \ 
/r(x2,  V2 )  J  \Xl  -  Xz  J 

0  \  (  yi  -  y2  \ 

fr{x 3,  J/3)  /  V^2  ~Xl  /  ’ 


and 

/r(*i,  !/i)  =  fir  +  /2rZi  +  hrVi  at  element  i, 

fr{x2 1  y2)  =  /lr  +  hrX2  +  /3ry 2  at  element  j, 

fr{x 3,  y3)  =  /lr  +  /2r^3  4-  /3ry3  at  element  J. 

Then  (6.5)  becomes 

<f={J/}T{ni}  +  {G}T{I>i}  +  y. 

The  potential  energy  of  the  normal  spring  is 


(6.11) 


P 

2 


P 

2 


V}T{2>;}  +  {G}T{G>}  +  y) 

}{ff}T{.D,}  +  {DJ}T{G}{G}T{flJ}  +  2{D.}T{/f}{G}T{DJ} 


+  2  {C.}T{J}  +  2  {C,}T{C}  + 


(6.12) 


Thus, 


p{H}{H}T=  \{H2}  {H2}t  { H3}t ),  (6.13) 

p{H}{G}T  =  |  {H2}  )({G!i}:r  {G2}T  {G^}7'),  (6.14) 

\Wj 

/{Gi}\ 

p{G}{H}T=  \{G2}  ({Fi}t  {H2}t  {H3}t),  (6.15) 

P{G}{G}T  =  |  {G2}  j  ({Gi}t  {G2}t  {G3}t  ) ,  (6.16) 

\{G3}) 

are  the  [Ke]  matrices;  and 


-p(^)w  =  -p(f)(|||).  (e,7) 

-P  (j)  fG)  =  “P  (f  )  ({§})’  (6'18) 

are  the  load  matrices. 

Then 


p{Hr}{Hs}T 

-> 

[^z(r)  *(*)]? 

r,s=  1,2,3, 

p{Hr}{Gs}T 

[A^(r)  j(5)], 

r,s  =  1,2,3, 

P{Gr}{Hs}T 

[A  j(r)  2(4)]) 

r,s  —  1,2,3, 

p{Gr}{Ga}T 

[7\;(r)y(5)], 

r,s  —  1,2,3, 

-"(f){*> 

{Fi(r)}> 

r  =  1,2,3, 

-"(f) 

— y 

{Aj(r) } , 

r  =  1,2,3, 
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where 


for  element  i 


i(l)  =  iu 
*(  2)  =  *2, 
*(3)  =  *3, 


and 


for  element  j 


j(l)  =  Ju 
j(  2)  =  j2, 

m=h- 


6.3  Shear  Contact  Matrix  for  Finite  Element  Covers 

In  Figure  3.8,  (xo,  yo)  is  on  the  edge  P2P3  and  is  the  assumed  contact  point  of  vertex 
Pi.  The  shear  spring  is  on  the  direction  P2P3,  and  connects  vertices  Pi  and  P0.  Let 


1  =  2  -  Z3)2  +  (y 2  -  ?/3)2, 


the  shear  displacement  of  P0  and  Pi  along  line  P2P3  is 

m 


(6.19) 


—  y  ( (zi  +  «i)  -  (x0  +  u0 )  (j/i  +  wi)  -  (yo  +  ^o) ) 


(x3  +  u3)  -  (x2  +  u2)  \ 
(2/3  +  V3)  —  (y2  +  v2)  )  ' 


Denote 


*S'o  =  (a:i-xo  yi-yo)(X 3  *2V 

\  2/3  —  2/2  / 

since  PqPi  is  small, 


(6.20) 


«  y  +  y  (*3  -x2  2/3  —  s/2 )  ^  ^ 

+  j(x2-X3  J/2  -  J/3  )  f  • 

I  \v0  / 


(6.21) 


W=  (kl)  =7Pi(*,’*,),r  (ll-ll)’  (6'22) 

{G>= (kl) =^k(io,!'o)it(^:^3)’  (623) 
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Denote 


where 


{Hr}  =  T 


{Gr} 


1  (  /r(xi,  yi)  0 

l  V  0  fr(xU  yi) 

1  f  fr(x 0,  yo)  0 

/  V  0  /r(^o,  yo) 


£3  -  Z2 
y3  -  y2 


X2  —  £3 


fr(x 0,  y0)  /  \  y2-y3 


and 

fr{xi ,  yi )  =  hr  +  hrX\  +  hry i  at  element, 
fr(xo,  y0)  =  /ir  +  /2r*o  +  /3r02  at  element j, 

Then  (6.21)  becomes 

d={H}T{D,}  +  {G}T{Di}  +  ^.  (6.24) 

The  potential  energy  of  the  shear  spring  is 

=  P^ 

2 

=  f  (W{A}  +  {G}rp,}  +  jj 

=  |  [{D.}7'{ff}{W)T{Di}  +  {D,}t{G}{G)t{D,}  +  2{0,}:r{ff}{G}:r{0,} 

x  2~i 

+  2  {fl.TW  +  2  (f )  {Dj}T{G}  +  (y)  J  '  (6'25) 


Thus, 

/{Hi}\ 

p{H}{H}T  =  \  {H2}  ({#i}T  {H2}T  {^a}T),  (6-26) 

\{H3}J 

f{Hi}\ 

p{H}{G}T=  {tf2}  ({Gi}r  {G2}t  {G3n,  (6.27) 

\{H3}J 

(  {Gi }  \ 

p{G}{H}T  =  {G2}  ({^i}T  {H2}T  {H3}T),  (6-28) 

f{Gi}\ 

p{G}{G}T=  \{G2}  ({G,}r  {G2}t  {G3}t),  (6.29) 

\{G3}J 


are  the  [Ke]  matrices,  and 


■'(?]<»> 


-(?)<«> 


(6.30) 


(6.31) 


are  the  load  matrices.  Then 


p{Hr}{Hs}T 

[Ki(r)  i(s)]> 

r,5  =  1,2,3, 

P{Hr}{Gs}T 

[Aj(r)  j(a)], 

r,  5  =  1, 2, 3, 

p{Gr}{Hs}T 

[A  j(r) 

r,5  =  1,2,3, 

p{Gr}{Gs}T 

r,s  =  1,2,3, 

-*(■  t)<*> 

-> 

}  1 

CO 

c^r 

1 — 1 

11 

-p  (j) 

{  ( r)  }  5 

r  =  1,2,3, 

where 

{  i(  1) 

=  *1, 

for  element  i 

i 

*(2) 

=  *2, 

1 

l  *(3) 

=  *3, 

and 

fi(i) 

=  ii, 

for  element  j 

< 

J(  2) 

=  J2, 

li(3) 

=  is- 

6.4  Friction  Force  Matrix  for  Finite  Element  Covers 

When  the  Coulomb’s  law  allows  sliding  between  two  sides  of  the  boundary  contacts, 
there  exists  the  friction  forces  in  the  two  sliding  sides  if  the  friction  angle  4>  is  not  zero. 

P\  is  on  element  i  and  P2,  P3,  Pq  are  on  element  j.  The  friction  force  is  calculated 
from  the  normal  contact  compressive  force  and  the  direction  of  the  friction  force  is  depend¬ 
ing  on  the  movement  of  Pi  relative  to  Pq  in  the  direction  from  P2  to  P3 .  Let  p  be  the  stiffness 
of  the  normal  contact  spring,  then  the  friction  force 
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T  =  p  ■  d  ■  s  ■  tan(<£), 


(6.32) 


where 

cl  is  the  normal  penetration  distance, 
tan (cf>)  is  the  friction  coefficient, 

s  -  sgn  (movement  of  Pi  relative  to  P0  in  the  direction  from  P2  to  P3) ,  and 


(  1,  if  x  >  0, 
sgn(x)  =  <  0,  if  x  =  0, 
l  —1,  if  x  <  0. 

The  friction  force  T  is  along  the  direction 

y((z3  -x2)  (y3  -y2)), 

and 

l  =  a/(x2  -  x3)2  +  (y2  -  y3)2 

is  the  length  of  line  P2P3. 

Then  the  potential  energy  of  friction  force  T  at  Pi  on  element  %  is 


„  7 , 

n/  =  t  (  Ul 


vi ) 


/  x3  -  x2  N 
\y3-y2  ) 


=  P{Di}T[Ti(x1,y1)]T 


(x3  -  x2)/l  \ 

(ys  -y2)/l ) 


(6.33) 


where 


{H} 


-l[Ti(xl,y1)]T 


f  x3  -x2\ 
\yz-V2 )  ' 


Then 


/{Hi} 

-P{H}  =  -7  {H2} 


(6.34) 


is  the  loading  matrix  for  element  i;  and 


~HHr} 


{Fi(r)  }>  r  —  1)  2,  3, 


where 


0  'I 

fr(x 1  yi) J 


fx3-x2\ 

\  2/3  —  s/2  y  * 


and 


fr{x  1  J/l)  —  /lr  +  f2rX l  +  fsrVl  • 


Similarly,  the  potential  energy  of  friction  force  T  at  P0  on  element  j  is 


n  /  = 


Vo  ) 


(  (*3  —  ar2)  \ 

\{yi-yi) ) 


Hd,}t{G], 


(6.35) 


where 


{G}  = 


][Tj(x0,y0)]T 


f  xz  -x2\ 

\y3-y2  J 


Then 

+T{G)  -  +T  I  {G2}  J 
\{Gz}J 


(6.36) 


is  the  loading  matrix  for  element  j;  and 


+JP{Gr}  ->  {-F?(r)}5  r  =  1,2,3, 


where 


-i  _  2.  (  fr{x 0  yo)  0  \ 

1  r)  “  /  V  0  fr{x 0  yo)  J 


(  X3  -  X2  \ 

\y3-y2 )  ’ 
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and 


fr{x o  yo)  =  fir  +  flr^O  +  fsrVO- 


6.5  Spring  stiffness 

The  stiffness  p  of  the  contact  springs  is  important.  If  p  is  too  small  the  penetration 
distances  become  too  large  so  that 

[1]  the  closed  contact  can  not  be  transferred  to  next  step; 

[2]  the  stresses  in  materials  can  be  reduced; 

[3]  deformation  along  joints  and  boundaries  can  be  wrong. 

If  spring  stiffness  p  is  too  large,  the  simultaneous  equation  can  be  nearly  linear  de¬ 
pendent  or  ill  conditioning  so  that 

[1]  the  solution  errors  can  be  unacceptable; 

[2]  the  iteration  method  may  not  converge; 

[3]  contact  displacement  may  not  be  correct. 

To  estimate  p,  the  relationship  of  spring  stiffness  p  and  material  Young’s  modulus 
E  is  computed.  Consider  a  rectangular  block  supported  by  two  springs  as  shown  by  Figure 
6.1, 
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the  block  stress  is 


2  F 
a  =  — 
a 

_a  _  2F 
E  aE 


db 


db 


2  bp 
aE 


(6.36) 


where 

a  is  the  length  of  the  rectangle 
b  is  the  height  of  the  rectangle 
F  is  the  value  of  two  vertical  point  loads 
a  is  the  vertical  stress  of  the  rectangle 
e  is  the  vertical  strain  of  the  rectangle 
db  is  the  vertical  displacement  on  the  top  of  rectangle 
d3  is  the  vertical  penetration  of  the  spring 
Assuming  a  =  b,  from  (6.36), 

db  =  ds (6.37) 

It  would  be  reasonable  if 

p  =  20  E  to  p  =  IOOjE 

is  chosen.  From  equation  (6.37),  the  block  displacement  db  is  40  times  to  200  times  of  spring 
displacement  ds.  Therefore  the  spring  displacements  are  neglected. 

In  the  current  program,  p  =  40 E  is  chosen.  There  are  several  ways  of  choosing  p 
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[1]  directly  enter  the  p  value; 

[2]  define  p  in  terms  of  E,  such  as  let  p  =  40 E; 

[3]  define  p  by  the  maximum  spring  displacement. 

6.6  Non-penetration  Physical  Springs 

The  options  of  non-penetration  springs  are  available  only  for  the  latest  discontinuous 
deformation  analysis  code.  The  algorithm  of  non-penetration  spring  is  still  in  testing  stage. 
This  algorithm  hopefully  will  be  transferred  to  the  manifold  method. 

For  the  real  materials,  the  contact  stiffness  is  not  as  large  as  the  program  had  set, 
but  the  penetration  still  keeps  small.  The  algorithm  of  the  low  stiffness  springs  with  small 
penetrations  will  make  the  contact  more  real  or  more  physical. 

The  previous  springs  are  mathematical;  while  the  new  springs  are  physical.  The  ad¬ 
vantages  of  the  physical  springs  are  in  the  following: 

fl]  the  spring  penetrations  are  small  so  that  the  contact  can  be  transferred  to  next  step 
without  missing; 

[2]  the  springs  can  have  the  physical  contact  stiffness; 

[3]  contact  springs  with  low  stiffness  will  form  better-conditioned  simultaneous  equa¬ 
tions; 

[4]  low  spring  stiffness  leads  faster  converge  of  iterative  solution; 

[5]  no  extra  unknowns  are  added; 

[6]  submatrix  structure  of  the  equations  are  reserved. 
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Chapter  7 

Two  Dimensional  Simplex  Integrations 


7.1  Simplex  Integration  on  a  Simplex 

For  the  manifold  method,  the  manifold  elements  are  the  integration  zones  which 
have  general  shapes,  therefore  the  integrations  are  more  difficult  than  the  integrations  of 
FEM.  However  analytical  solutions  were  found  for  many  cases  of  manifold  method.  The 
finite  element  method  computes  the  integrations  of  complex  functions  in  simple  domains; 
the  manifold  method  computes  the  integrations  of  simple  functions  in  complex  domains. 
The  integrations  on  complex  domain  can  be  reduced  to  the  integration  on  the  simplex.  The 
simplex  is  the  oriented  simplest  domain. 

The  simplex  has  the  most  simple  shape  in  1, 2, 3, . . . ,  n  dimensional  space.  Differ¬ 
ent  from  the  ordinary  integration,  the  simplex  integration  has  only  the  simplex  as  the  integral 
domain.  The  simplex  also  has  positive  or  negative  orientations.  Positive  or  negative  orien¬ 
tations  define  positive  or  negative  volumes  respectively.  Figure  7.1  shows  the  0,  1,  2,  3 
dimensional  simplex. 


Figure  7.1 .  0,1 ,2,3  dimensional  implex 
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The  0  dimensional  simplex  is  a  point  P0,  its  volume  is 


—  |  11  =  1 
0!  1  ' 


(7.1) 


The  0-d  simplex  integration  can  be  considered  as  a  normal  real  number. 
The  one  dimensional  simplex  PqP\  is  an  oriented  segment,  its  volume  is 


1 

V. 


1  Xio 

1  xu 


^11  —  ^10 


(7.2) 


The  volume  of  the  simplex  P\Pq  is  the  negative  volume  of  the  simplex  PqP\.  The 
ordinary  1-d  integrations  are  the  simplex  integrations: 


Therefore,  for  any  co-line  point  P2, 


r  P2 
'Po 


fPi  fPi  rP2  fPi  rPi 

f(x)dx+  /  f(x)dx  =  /  f(x)dx+  /  f(x)dx+  /  f{x)dx  =  /  f(x)d 
JPi  JPo  JPi  Jp2  JPo 


The  1-d  integration  addition  is  the  same  as  vector  addition.  The  integrations  on  neg¬ 
ative  vectors  and  positive  vectors  can  be  nullified. 


Two  dimensional  simplex  P0P1P2  is  an  oriented  triangle,  its  volume  is 


1 

2! 


1  Xio  £20 
1  xu  £2i 
1  Z12  Z22 


(7.3) 


The  volume  of  simplex  PiPqP2  is  the  negative  volume  of  simplex  PqPiP2. 


7.2  Simplex  Integration  on  a  General  Polygon 

Unfortunately,  the  two  dimensional  ordinary  integrations  are  volume  integration, 
where  the  volume  is  always  positive.  For  the  2  dimensional  ordinary  integration,  the  inte¬ 
gration  domains  have  no  orientations,  therefore  the  integrations  have  no  algebraic  addition 
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on  oriented  domains.  It  is  then  necessary  to  introduce  simplex  integration  on  a  simplex, 
which  is  a  simple  oriented  domain. 

The  simplex  integration  is  not  intended  to  do  integrations  on  a  simplex  only.  Obvi¬ 
ously  a  complex  shape  always  can  be  subdivided  into  simplex,  so  simplex  integration  can  be 
computed  in  each  simplex,  the  summation  of  simplex  integrations  is  the  ordinary  integration 
over  the  complex  shape.  However  this  is  also  not  the  way  of  using  simplex  integration. 

Giving  a  polygon  PiP2P2PiP$P(i  with  P6  =  Pu  such  that  P,  rotates  at  the  same 
direction  as  from  ox  to  oy.  For  any  point  P0,  the  algebraic  addition  of  the  2-d  simplex  vol¬ 
ume  (area)  P0P1P2,  P0P2P3,  P0P3P4,  P0P4P5  and  P0P5P1  is  the  area  A  of  the  polygon. 
Let  P0  =  (0, 0), 


5 

1 

0 

0 

1  5 
=  5^ 

E 

1 

1 

X\  j 

£2  i 

X\  1  *^2  i 

^1  i+1  i+1 

t=l 

I 

*^1  t+ 1 

i+1 

1  =  1 

(7.4) 


Figure  7.2  shows  the  addition  of  plus  and  minus  areas  of  simplex.  In  Figure  7.2,  the 
area  of  simplex  P0P2P3,  P0P3P4,  P0P4P5  and  P0P5P1  are  positive;  the  area  of  simplex 
P0P1 P2  is  negative.  The  algebraic  sum  is  exactly  the  area  of  polygon  Pi  P2 P2 P4  P5 P6 .  The 
area  A  is  then  represented  by  the  coordinates  of  boundary  vertices. 


Figure  7.2.  Addition  of  plus  or  minus  area  of  simplex 
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In  general,  the  simplex  integration  can  compute  ordinary  integrations  without  subdi¬ 
viding  2-d  domains  to  triangles.  FEM  mesh  is  unnecessary  for  simplex  integration.  Using 
simplex  integration,  the  integration  of  any  polynomials  can  be  represented  by  the  coordi¬ 
nates  of  boundary  vertices  of  generally  shaped  polyhedron. 


7.3  Two  Dimensional  Simplex  Integrations 

The  two  dimensional  integrations  here  in  this  section  are  used  for  two  dimensional 
manifold  method  algorithm.  Since  the  displacement  function  is  linear  function  of  coordi¬ 
nates  (x,  y),  the  integrands  are  of  degree  0,  1,2. 

A  2  dimensional  simplex  has  3  vertices 

Po,P\ i  P7- 


P0  :  (xio,  £20) 

Pi  ■  (£11,  £21)  (7-5) 

P2  :  (£12,  £22) 

The  2  dimensional  coordinate  simplex  has  3  vertices 

U0,UUU2. 


Uq  :  (  0,  0) 

Ui  :  (  1,  0)  (7.6) 

U2:  (  0,  1) 

The  following  coordinate  transformation 


{ui,u2)  (xi,x2) 


transfers  coordinate  simplex 


U0U1U2 
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to  normal  simplex 


P0P1P2 


X1  —  ZlO  (1  — ^lu*)+  XH  U1  +  x12  U2 

x2  =  x20  (1  ~  Ui)+  X21  U\  +  #22  U2 

The  Jacobi  determinant  is 

J  _  D(Xh  X2) 

D(u  1,  U2 ) 

dx\  dx\ 

_  dui  du2 

dx2  dx2 
du 1  du2 

__  ^12-^10 
#21  ~~  £20  ^22  -  X20 

1  #10  #20 
=  1  xU  ®21 

1  X12  X22 


(7.7) 


(7.8) 


Since  P0PiP2  is  a  2-dimensional  simplex  with  non-zero  volume,  J  is  non-zero. 
Translation  can  be  rewritten  as 


(uo,Ui,U2)  -4  (1, #1,0:2) 


1 

= 

u0 

+ 

U\ 

+ 

U2 

X\ 

= 

Uo 

+ 

#ii 

Ul 

+ 

X12 

V>2 

*2 

=  ^20 

Uq 

+ 

x2\ 

U\ 

+ 

X22 

Two  dimensional  simplex  integration 

Sp0P1P2{rni,m2) 


on  a  2  dimensional  simplex 


P0P1P2 


is  defined  as  normal  integration  times  the  sign  of  determinant  J. 
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Sp0p1Pi(mi,m2) 

=  sign(J)  f  I  x™1  x™2  dx\dx2 

J  JP0P1P2 

=  sign(J)  f  I  x’P'x™2 \ J\du\du2 

J  JU0U1U2 

=  J  j  j  x™lx™2du\du2 

J  JU0U1U2 


x™1  x™2duidii2 


(7.10) 


XT‘  ~  (!>**») 


then  the  two  dimensional  simplex  integration  can  be  represented  by  the  following  basic 
forms  of  simplex  integration: 


S\  =  J  J  U%q  u'l  u'f  du\dll2 

J  Ju0ulu2 

=  I  UlQu\1Ul2duidu2 

J  Juo,Ui  ,«2>0 

nl“«2 

Ul2  u'l  UlQ  dxi\du,2 

_ 

=  J  u\ 2  (^J  uj*(l  —  u2  —  r£i)l0dui^  du2  ( 

After  ii  times  of  integration  by  parts,  the  inner  integration  can  be  computed. 

rl~u  2 

/  Uj1  (1  —  U2  —  UiY°dui 

Jo 

rl-u2  I 

=  /  di—ru^l-ut-u.r 

Jo  *i  +  l 


(7.11) 


^l+l 


*1+1/1  \  *0  *  u2 
u 1  (1  -  u2  -  Ui)  0 


I-1-U2  I 

-  U;>+V((l  -  u2  -  Ul)«) 

Jo  *1  +  -t 

r  1— »a  -j 

=  -/  -«!)-) 

Jo  *1+1 
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1  —  U2 


^0  21  +  1  /  1  \  t*0  “1 

~  .  7«i  (1  “  U2  -  Ui) 

*i  +  l 


C?Ul 


*o(*0  —  1) 


\io-2 


0 

'1  —  «2 

y.r»  i  v.rk  —  I  I  •  i  o 

r  .  i \r  .  0\ui1+2(1  ~  u2  -  «i)‘° 
/o  (*i  +  l)(*i  +  2) 

’1_“2  jo(*o  -  l)(io  -  2)  „ll+3, 


/' 

f 


uj1+i5(l  —  U2  ~  Ui)*°  3dUi 
_ ,,*l+*0-l 


(*i  +  l)(*i  +  2)(*1  +  3) 

1-“2  *o(*o  —  1)(*0  —  2) ...  2 


(*i  +  l)(*i  +  2)(ii  +  3) . . .  (t"i  +  io  —  1)  1 


(1  —  U2  —  u\)1  du\ 


f  1-U2 

Jo 


*o(*o  —  1)(*0  —  2) ...  1 


(*i  +  l)(*i  +  2)(*i  +  3) . . .  (*i  +  *o) 

/*1-“2  i  .  . 

/  — — P * 

Jo  (ti+*o)!  1 


u[1+iodui 


» 1  —  ti2 


*o!*i! 


(*i  +  *o  +  1)! 


(/(Uj 


d  “Ho+1  > 


(*o+*i  +  l)!Ul 


*0+*l  +  l  |1-“2 

lo 


(ii  +  io  +  1)! 


(7.12) 


u^du 


1  — «2 

2  /  Uj1  (1  —  U2  —  U\)l°du 


1 


/ 


'0 


-  °2) 


to!ti!t*2! 


io  +  il  +  Vuo 


l0!*i! 


(*o  +  *1  +  1)! 


(*o  +  *1  +  *2  +  2)! 


(7.13) 


Based  on  the  previous  formula  of  S i,  the  two  dimensional  integration  of  polynomi¬ 
als  of  degree  0,  1,  2  can  be  computed. 


SVoPiP2(0>0)  =  sign(J)  /  /  dx1dx2 

J  Jp0PiP2 

—  J  I  / 

J  JUoUxlh 

SPoPiPvi®^)  =  sign(J)  /  /  dx\dx2 

J  JP0P1P2 
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(7.14) 


=  J- 


0! 


=  Ij 


(0  +  2)!  2 

SpoPiPii1’®)  =  sign(J)  /  /  xidxidx 2 

•/  d  Pq  P !  j 


L  P'2 


-'ll 


UoUtUi 


(x10u0  +  X11U1  +  x12u2)du1du2 


J- 


1! 


r(*io  +  *n  +  x12) 


(1  +  2)! 

Spop1p2(l,0)  =  sign(J)  /  x\dx\dx2 

J  J  Po  Pi  P 2 


1 

=  qJ(xio  +  Xn  +Xu) 


Sp0Pi  p2(0, 1)  =  sign(J)  /  /  x2dxidx2 

J  J  Pq  Pi  P% 


—  qJ(X  20  +  2-21  +  ^22) 


^oPiPa(2>0)  =  sign(J)  /  /  x\dx-ydx2 

J  J  Pc  Pi  Pi 

-'ll 


’U0UtU2 


Po  Pi  P-2 

(x10u0  +  xnui  +  X\2u2)2  du\du2 


=  J 


+  J 


2! 


(2  +  2)! 

1! 


(*10 +  *ii  +*12) 


(2  +  2)! 

(0  +X10X11  +XjoX^2 

+X11X10  +0  +X11X12 

+X12X10  +X12X11  +0) 

Sp0Pi  p2(2,  0)  =  sign( J)  /  /  x\dx^dx 2 

J  JPoPlPl 


J 


1_ 

24 


(2xioXio  +Xio^n  +^10^12 
+X11X10  +2xuxn  +X11X12 
+  X12X1Q  ~\~x12xll  +2^12X12) 


Spop1p2(0^  2)  =  sign{  J)  /  /  x\dx\dx2 

J  JP0P1P2 


4 


(2X20^20  +  X2oX2i  +X20X22 

+  X21X20  +2x2lX21  +X21X22 

+X22X20  +X22X21  +2X22X22) 


(7.15) 


(7.16) 


104 


SpQ Pi p2 ( 1  >  1)  —  sign^J')  /  j  X\X2dX\d,X 

J  JP0P1P2 

-'JL 


Po  Pi  ^2 

{xioUo  +  Xu  «1  +  ^12^2) 

U0UlU2 

(X20U0  +  a;2i«i  +  X2iuz)du\du2 
2! 

(Z10Z20  +  X11S21  +  £12^22) 


=  J 
+  J 


(2  +  2)! 

1! 


(2  +  2)! 

(0  +Z10Z21  +Z10Z22 

+  ^11^20  +0  -\-XuX22 

+^12^20  +^12^21  +0) 


Sp0Pi  P2(l,l)  =  sign(J)  /  /  xixzdxidxi 

J  Jp0p1p2 


=  —  J(2xiox2o  +  2xnx2i  +  2x12X22 
+  £10^21  +  210^22  +  a;  11  ^20  +  3:11X22  +  £12^20  +  3:122:21) 


(2xioX20 

+3:11X20 

+X12X20 


+X10X21  +3:103:22 

+2xhx2i  +X11X22 
+X12X21  +2X12X22) 


(7.17) 


7.4  Simplex  Integration  for  Two  Dimensional  Manifold  Method 

Since  simplex  integrations  always  have  the  Jacobian  J  as  factor  and  J  is  an  oriented 
area,  the  integrations  on  positive  area  and  negative  area  can  be  neutralized.  Denote 


PlP2iV--Pn  Pn+l=Pl 


with 


P;  =  (xij,X2f) 


as  a  polygon,  the  integrations  (7. 14)-(7. 17)  on 

Pi  p2  Pi  ‘  '  '  Pn 
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are  computed  by  summations.  (7. 1 8)-(7.2 1)  are  the  integrations  of  manifold  method. 
Here  P0  can  be  any  point,  Po  =  (0, 0)  are  chosen  in  order  to  have  simpler  formulae.  Com¬ 
puted  by  simplex  integrations,  integrals  (7.18)-(7.21)  are  represented  by  the  coordinates  of 
the  boundary  vertices  only. 


:1  ^PoP^Pm-iCM) 

—  1  Xl  k  x 2 

-  2  2s k= i  X\  fc+i  x2  I 

/  /  xldxidx2  =  V\  ^PoPfcPfc+i(l50) 

J  J (A)  *-1 

[  [  x2dxidx2  =  7  V'  , 

7  7(a)  6  z"^fc=1 

/  J  xldxidx2  =Yll=  _ j  '5p0P*P*+i(2, 0) 

1  \ 

=  12  ^*=i 

(x?  it  +«i  *+i  +*i  Jtxi  *+i) 

x]dx\dx2  —  — 

1  24  ^fc=l 


x\  k  x2  k  ,  ,  \ 

(#1  Jb  +  ^1  fc-fl  J 

fc-f-l  ^2  k -4*1 

X1  k  x2  k  f  .  \ 

(^2  k  +  x2  k+ 1  j 

A: H- 1  *^2  &-fl 


k  x2  k 
x  1  A-f-1  *^2  fc-fl 


(7.18) 


(7.19) 


[  [  xidxidx2  —  —  V'' 
J  J  (A)  2  24 


Xx  X2 

Xl  fe+1  X2  k+1 

(2xi  kX  1  fc  +Xi  fcXi  fc+1 

+xx  fe+ixj  fc  +2xi  fc+ixj  t+x) 
XI  k  X2  k 

Xl  k+1  X2  k+1 

(2x2  kX2  k  +X2  kX 2  fc-fl 

+X2  it+1^2  k  +2X2  k+  lx2  fc+l) 


(7.20) 


/  /  xix2dx\dx2  =  „  5P0PfcPfc+i(l>  1) 

7  7(A)  ^k=i 


1  v — > n 

24  £-^k—  1  J  x  rv  7  j.  n>  (  J-  , 

(2xi  fcx2  it  +  2xi  jt+ix2  *+i  +  xx  fcx2  jfc+i  +  xx  fc+ix2  k) 

I  I 


X1  k  x2  k 
X1  x2  Ar-j-l 


(A) 


1  V — "V  n 

xix2dx1dx2  =  —  >  , 

24  *-^k=i 


Xx  X:  X2 

Xl  it+1  X2  x+l 

(2xx  fcx2  x  +Xx  fcX2  fc+i 

+xi  fc+xx2  k  +2xx  fc+ix2  t+i 


(7.21) 
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Chapter  8 

Equation  Solver  and  Open-Close  Iterations 
for  Inertia  Dominant  Equilibrium  Equations 

8.1  Time  Step  Based  Large  Displacement  Analysis 

The  manifold  method  computations  follow  the  time  steps.  The  time  steps  are  small 
enough  that  the  second  order  displacements  are  neglected.  All  the  geometric  and  physical 
parameters  have  to  be  transferred  from  the  end  of  the  previous  time  step  to  the  beginning  of 
the  next  time  step.  The  following  items  are  to  be  transferred: 

[1]  stresses  of  each  element, 

[2]  strains  of  each  element, 

[3]  velocities  of  each  element, 

[4]  geometry  of  the  joint  boundaries  and  elements, 

[5]  all  closed  contacts 

The  position  and  state  parameters  of  the  closed  contacts  have  to  go  to  the  next  step. 
The  geometric  parameters  and  physical  parameters  of  contacts  will  be  transferred: 

[1]  the  contact  vertex  and  edge, 

[2]  the  position  of  contact  point, 

[3]  the  normal  displacement  and  normal  force, 

[4]  the  shear  displacement  and  shear  force, 

[5]  locking  or  sliding  as  contact  state. 

The  manifold  method  computes  both  statics  and  dynamics  by  time  steps.  For  large 
displacements  and  deformations,  the  statics  is  the  ultimate  stabilized  state  of  dynamics  af- 
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ter  sufficient  long  times.  Such  a  stabilized  state  can  be  reach  only  if  some  kind  of  energy 
reduction  is  counted  in  the  computation. 

The  previously  customized  one-step  static  computation  is  correct  only  if  the  second 
order  displacements  can  be  omitted.  In  such  a  case,  physically  one  infinite  large  time  step  A 
is  chosen  under  the  small  ultimate  displacements.  The  inertia  force  would  be  proportional 
with 


The  value  of  (8.1)  is  zero  if  time  step  A  is  very  large,  then  the  inertia  force  is  zero. 
This  is  the  reason  why  the  one  step  statics  do  not  consider  the  inertia  force  or  mass  matrices. 

The  current  version  of  the  manifold  code  is  a  simple  version.  In  the  beginning  of  a 
time  step,  the  dynamic  computations  inherit  the  velocity  of  the  end  of  the  last  step.  For  static 
computation  it  is  assumed  that,  the  initial  velocity  of  each  time  step  is  zero.  The  current 
static  computation  is  obviously  incomplete,  the  way  of  energy  reduction  has  to  be  studied 
in  the  future. 


8.2  Open-Close  Iterations 

Within  each  time  step,  the  global  equations  have  to  be  solved  repeatedly  while  se¬ 
lecting  the  lock  positions.  The  procedure  of  adding  and  removing  stiff  springs  is  open-close 
iteration. 


If  a  contact  has  a  tensile  contact  force  from  the  normal  spring,  the  two  sides  will 
separate  after  the  removal  of  this  stiff  spring.  If  the  vertex  penetrated  the  edge  in  other  side 
of  the  contact,  a  stiff  spring  is  applied. 


For  each  contact,  there  are  three  modes: 
mode  changing  is  in  the  following, 
mode  change  condition 

open  -  open  N  >  0 

open  -  slide  N  <  0  and  |T|  > 

open  -  lock  N  <  0  and  |T|  < 

slide  -  open  N  >  0 


open,  sliding  and  lock.  The  criteria  of  the 


tancf>\N\ 

taruf)\N\ 
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slide  -  slide 

N  <  0  and  1 1|  —  / 

slide  -  lock 

N  <  0  and  t  \\f 

lock  -  open 

N  >  0 

lock  -  slide 

N  <  0  and  |T|  >  tan4>\N\ 

lock  -  lock 

N  <  0  and  |T|  <  tan(f)\N\ 

where  p  is  spring  stiffness 

n  is  normal  displacement  vector  pointing  vertex 

N  is  normal  displacement,  N  >  0  is  open 

t  is  shear  displacement  vector  pointing  P2-P3 

T  is  shear  displacement,  T  >  0  if  t  in  same  direction  as  P2P3 

||  means  two  vectors  point  the  same  direction 

t  is  friction  force  vector  pointing  P2P3 

(j)  friction  angle 

Finding  the  mode  changing,  the  following  operations  will  be  done: 


mode  change 

operation 

open  -  open 

no  changing 

open  -  slide 

apply  pair  of  friction  forces 

open  -  lock 

apply  normal  and  shear  springs 

slide  -  open 

delete  friction  forces 

slide  -  slide 

no  changing 

slide  -  lock 

delete  friction  forces  and  apply  shear  spring 

lock  -  open 

delete  normal  and  shear  springs 

lock  -  slide 

delete  shear  spring  and  apply  friction  forces 

lock  -  lock 

no  changing 

The  open-close  iterations  to  ensue 

[1]  no-penetrations  in  the  open  contacts, 

[2]  no-tensions  in  the  contacts  with  normal  springs. 
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The  two  conditions  have  to  be  fulfilled  in  all  contacts.  If  the  two  conditions  are  not 
fulfilled  after  five  times  of  open-close  iterations,  the  time  step  will  be  reduced  to  one  third, 
and  the  open-close  iteration  continues. 


8.3  SOR  Iteration  Method 

The  abbreviation  “SOR”  stands  for  successive  over  relaxation  method.  This  method 
is  for  solving  linear  equations, 


=  [F] 


(8.2) 


Equation  (8.2)  has  submatrix  structure 


An 

A12 

A13 

A14 

•  *  ■  A\ n  ^ 

( *i\ 

(F1\ 

A21 

A22 

A23 

A24 

A2  n 

f2 

A31 

A32 

A33 

A34 

Ain 

X3 

Ft 

A41 

A42 

A43 

A44 

.  .  .  A^n 

x4 

Fa 

\  AjjI 

Aji2 

An3 

Ajj4 

•  •  •  Ann  / 

\xn) 

\FnJ 

(8.3) 


where 

[1]  [A]  is  a  n  x  n  coefficient  matrix, 

[2]  [X]  is  a  n  x  1  unknown  matrix  , 

[3]  [F]  is  an  x  1  matrix  of  free  terms.  The  elements  of  these  matrices  are  still  subma¬ 
trices. 

[1]  The  elements  Atj  of  matrix  [A]  are  q  x  q  matrices, 

[2]  The  matrix  [A]  is  symmetric,  [Aij}T  =  [A,-,-] 

[3]  The  elements  Xt  of  matrix  [X]  are  q  x  1  matrices, 

[4]  The  elements  F;  of  matrix  [F]  are  q  x  1  matrices. 
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Denote  the  diagonal  matrix  [D]  as 


Then  there  are 

«UJ  -  PI)  +  p])M  =  p] 

pim  =  pi  -  m  -  pum 

[X]  =  pj-'p]  _  p]-*([A]  -  [£])[*] 

m  =  [Gi-ppi  (s.4) 

where 

[G\  =  WW) 

[B)  =  [D)-1([A}-[D]) 

[1]  [B]  is  a  n  x  n  coefficient  matrix, 

[2]  [G]  is  a  n  x  1  force  matrix , 

[3]  [X]  is  a  n  x  1  unknown  matrix ,  The  elements  of  these  matrices  are  still  submatrices. 

[1]  The  elements  Du  of  matrix  [D]  are  q  x  q  matrices, 

[2]  The  elements  B{j  of  matrix  [B]  are  q  x  q  matrices, 

[3]  The  elements  G,  of  matrix  [G]  are  q  x  1  matrices. 

[-Dji]  =  [.Ajj], 

Po)  =  |0],  l.fl  7=  j 

<  Pol  =  piil-'py]. 

Pyl  =  Piil  ‘[/M,  'S' 

l.  [B„]  =  [0], 


HI 


For  SOR  method  the  factor  of  over  relaxation  is  u>, 

1  <  u  <  2 

Denote  is  the  solution  after  m  iterations,  then  the  solution  of  the 

next  iteration  is 
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8.4  Simple  Iteration  Method 

The  simple  iteration  method  of  equation 

[A)[X]  =  [F\ 

has  a  simpler  algorithm.  Generally  speaking,  the  simple  iteration  method  is  mush  more 
slower  than  SOR  iteration  method.  The  formula  of  simple  iteration  method  is  mathemat¬ 
ically  clear  and  clean.  However,  practically  this  method  never  became  a  major  equation 
solver,  because  of  its  low  efficiency  comparing  with  the  other  iteration  methods.  It  seems 
this  method  will  stay  in  textbook  forever  to  give  the  beginners  a  fresh  idea  about  what  is 
iteration  method. 

It  is  often,  that  things  can  change  in  an  unexpected  manner.  Two  thing  happened 
recently: 

[1]  parallel  computation, 

[2]  inertia  dominant  computation  in  discontinuous  deformation  analysis. 

Now,  simple  iteration  method  is  the  equation  solver  which  is  most  suitable  to  the 
parallel  computation.  In  the  mean  time,  the  diagonal  dominant  matrices  of  DDA  and  mani¬ 
fold  method  are  good  enough  to  use  the  simple  iteration  method.  Choosing  small  time  steps, 
the  coefficient  matrices  can  be  almost  diagonal  matrices,  then,  there  will  be  almost  no  dif¬ 
ference  between  simple  iteration  method  and  SOR  iteration  method. 

As  before,  the  equation  is 


/An 

A 12 

-^13 

A14 

•  •  •  A\n  ^ 

/Xi\ 

(F 1\ 

A21 

A22 

A23 

A24 

A.2n 

x2 

f2 

^•31 

A32 

A33 

A34 

A3n 

X3 

F3 

A41 

A42 

A43 

-4.44 

•  •  ■  A.4n 

X4 

F4 

V  Anl 
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a' 
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\Xn) 

\fJ 

The  algorithm  of  the  simple  iteration  method  is 

/x\m+1)\ 

\r(m~ f  1) 

v(m+!) 

a3 


\X{nm+1)/ 


113 


8.5  Time  Step  Algorithm  and  Iteration  Convergence 

The  principles  of  choosing  time  steps  A  are 

[1]  A  is  small  enough,  so  that  the  second  order  displacements  are  neglected; 

[2]  A  is  small  enough,  so  that  the  SOR  iteration  will  converge  in  less  than  30  iterations; 

[3]  A  is  small  enough,  so  that  the  open-close  iterations  will  converge  in  less  than  6  it¬ 
erations, 

[4]  A  is  large  enough,  so  that  the  computation  will  represent  larger  time  span  and  the 
displacements  are  stabilized  if  possible. 

There  are  three  options  to  define  time  step: 

[1]  directly  enter  the  time  step, 

[2]  enter  the  allowable  maximum  step  displacement,  then  define  time  step  from  the  max¬ 
imum  allowable  step  displacement, 
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[3]  compute  the  maximum  allowable  step  displacement,  then  define  time  step  from  the 
maximum  allowable  step  displacement. 


Denote 

T>  as  the  maximum  allowed  step  displacement, 

U  as  the  mass  of  the  average  elements, 

T  as  the  maximum  load  on  the  average  elements, 
V  as  the  maximum  velocity  of  elements, 

The  equation  of  time  step  A  is 


let 


£A2  +  VA-£>  =  0 


A  = 


^(-V  +  \/V2+4  QV) 
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Chapter  9 

Applications  of  Numerical  Manifold  Method 


9.1  Joint  Computations 

Figure  9.1  shows  the  ability  of  the  manifold  method  to  compute  a  joint  or  fracture. 
Figure  9.2  shows  the  deformations  of  a  domain  with  a  self  intersected  curved  joint 

9.2  Block  Computations 

Figure  9.3  shows  the  failure  of  an  arch  under  the  point  load  on  the  center  block  and 
self  weight. 

Figure  9.4  shows  the  failure  of  a  gravity  dam  with  rock  foundation.  The  loads  are 
the  upstream  water  pressure  and  the  self  weight  of  the  dam. 


9.3  Slope  Sliding 

Figure  9.5  is  the  result  of  slope  sliding  of  rock  blocks.  It  can  be  noticed  that,  the 
center  block  separated  two  adjacent  blocks  during  the  sliding.  The  result  is  consistent  with 
the  lab  test. 

Figure  9.6  is  a  soil  slope  which  slides  along  a  circular  surface.  The  circle  sliding 
computation  satisfies  all  equilibrium  conditions. 


9.4  Failure  of  Structures 

Figure  9.7  shows  the  deformation  of  the  joints,  blocks  and  the  continuous  materials 
of  double  simply  supported  beams 

Figure  9.8  shows  the  deformation  of  the  joints,  blocks  and  the  continuous  materials 
of  jointed  double  cantilever  beams 
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FIGURE  9.1  Single  joint 


I 
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FIGURE  9.2  Deformations  of  a  domain  with  a  complex  joint 


FIGURE  9.3  The  failure  of  a  deformable  arch 
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FIGURE  9.4  Failure  of  dam  and  rock  foundations 
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FIGURE  9.5  Multi-block  sliding 
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FIGURE  9.6  Circular  sliding 
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FIGURE  9.7  Deformation  of  double  simply  supported  beams 


125 


FIGURE  9.8  Deformation  of  jointed  double  cantilever  beams 


9.5  Conclusions 


This  new  theory,  entitled  the  Manifold  Method  of  Material  Analysis,  combines 
physical  meshes  and  mathematical  meshes.  These  physical  meshes  provide  the  means  to 
consider  both  jointed  and  continuous  materials,  and  even  different  material  phases  (i.  e. 
solid,  gas  or  liquid).  At  present,  a  theory  for  the  manifold  method  has  been  accomplished, 
as  has  a  first  generation  2-D  dynamic  computer  code.  The  preliminary  results  are  encour¬ 
aging  ( for  example,  the  convergence  of  the  solutions  has  been  established).  Finite  element 
and  DDA  formulations  are  special  cases  of  this  developing  theory.  A  brief  listing  of  a  few 
of  the  advantages  of  the  manifold  method  follows: 

[1]  free  surface  and  flexible  boundaries 

[2]  analysis  not  hindered  by  boundary  conditions 

[3]  free  form  elements  (any  shape) 

[4]  conservation  of  energy 

[5]  obeys  Coulomb’s  law 

[6]  very  small  to  very  large  deformation 

[7]  statics  and  dynamics  possible 

[8]  analytically  correct 

[9]  continuous  and  discontinuous  analysis 
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Appendix:  The  mathematical  manifold 


The  “  Manifold  ”  is  a  main  subject  of  differential  geometry  and  topology  in  modern 
mathematics.  The  basic  structure  of  the  manifolds  is  a  finite  cover  system  and  the  connec¬ 
tion  functions  among  the  covers.  Different  from  the  mesh,  the  covers  can  be  overlapped  or 
folded  to  present  a  combinatorial  space.  Therefore  the  manifold  method  is  to  pursue  global 
solutions  on  global  spaces.  There  were  no  direct  connections  before  between  manifolds  and 
engineering  analysis.  The  essential  difficulty  is  the  complicated  boundaries  and  discontin¬ 
uous  interfaces  of  the  real  engineering  cases. 

Under  the  generalized  definition  of  manifold,  there  are  two  independent  mesh  sys¬ 
tems  now:  mathematical  mesh  and  physical  mesh.  The  mathematical  mesh  consists  of  the 
folded  covers.  Independent  cover  functions  are  defined  on  each  cover.  The  weight  functions 
connect  all  locally  defined  functions  to  a  combinatorial  global  function. 
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Simplex  Integration  for  Manifold  Method 
and  Discontinuous  Deformation  Analysis 


Gen-hua  Shi 

Geotechnical  Lab,  US  Army  Engineer  Waterways  Experiment  Station 
Vicksburg,  MS  39180-6199 

Abstract 

At  least  every  engineer  has  to  compute  the  volume  of  generally  shaped  blocks.  Is 
there  a  formula  where  the  volume  is  precisely  represented  by  the  coordinates  of  boundary 
vertices?  If  block  movements  are  considered,  the  center  of  gravity  has  to  be  computed?  Is 
there  a  formula  where  the  center  of  gravity  is  also  represented  by  the  coordinates  of  bound¬ 
ary  vertices?  The  simplex  integration  developed  for  DDA  computation  can  also  solve  these 
questions.  The  convergency  and  accuracy  of  DDA  algorithms  depend  upon  mainly  the  an¬ 
alytical  integrations  on  complex  shapes.  Simplex  integrations  are  accurate  solution  on  n- 
dimensional  generally  shaped  domains.  The  integrand  could  be  any  n-dimensional  polyno¬ 
mials.  The  DDA  computations  of  three  rock  failure  cases  are  presented. 


Simplex  Integration  on  a  Simplex 

Simplex  has  the  most  simple  shape  in  1, 2, 3, . . . ,  n  dimensional  space.  Different 
from  the  ordinary  integration,  simplex  integration  has  only  simplex  as  integral  domain. 
Simplex  also  has  positive  or  negative  orientations.  Positive  or  negative  orientations  define 
positive  or  negative  volumes  respectively. 


Pa 


Figure  1 .  0,1 ,2,3  dimensional  implex 
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The  0  dimensional  simplex  is  a  point  Po,  its  volume  is 


1 

0! 


1 1  =  1 


The  0-d  simplex  integration  can  be  considered  as  a  normal  real  number. 


One  dimensional  simplex  PqP\  is  an  oriented  segment,  its  volume  is 


1  £io 
l!  1  £n 


Xu  -  £l0 


The  volume  of  the  simplex  PiPo  is  the  negative  volume  of  the  simplex  P0Pi  .  Ordinary  1-d 
integrations  are  simplex  integrations: 


f{x)dx 


Therefore,  for  any  co-line  point  P2, 


r  P2  r  Pi  f  Pi  f  Pi  r  Pi 

/  f(x)dx-j-  /  f(x)dx  =  /  f(x)dx- b  /  f(x)dx-\-  /  f(x)dx 

Jpo  J  P2  Jpo  JPi  'Pi 


The  1-d  integration  addition  is  the  same  as  vector  addition.  The  integrations  on  negative 
vectors  and  positive  vectors  can  be  nullified. 


Two  dimensional  simplex  P0PiP2  is  a  oriented  triangle,  its  volume  is 

1  X\Q  £20 

1  xn  x2i 
1  £12  £22 

The  volume  of  simplex  P1P0P2  is  the  negative  volume  of  simplex  P0P1P2. 


1 

2! 


Three  dimensional  simplex  P0P1P2P3  is  a  oriented  tetrahedron,  its  volume  is 


1 

XlO 

£20 

£30 

1 

1 

Xu 

X21 

£31 

3! 

1 

X\2 

X22 

£32 

1 

£13 

^23 

£33 

The  volume  of  simplex  P1P0P2P3  is  the  negative  volume  of  simplex  P0P1P2P3. 

Unfortunately,  two  dimensional  and  three  dimensional  ordinary  integrations  are  vol¬ 
ume  integration,  where  the  volume  is  always  positive.  For  the  2  or  3  dimensional  ordinary 
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integration,  the  integration  domain  has  no  orientation,  therefore  the  integration  has  no  alge¬ 
braic  addition  of  oriented  domain.  It  is  then  necessary  to  introduce  simplex  integration  on 
a  simplex,  which  is  a  simple  oriented  domain. 

The  simplex  integration  is  not  intended  to  do  integrations  on  simplex  only.  Obvi¬ 
ously  a  complex  shape  always  can  be  subdivided  into  simplex,  so  simplex  integration  can 
be  computed  in  each  simplex,  the  summation  of  simplex  integrations  is  the  ordinary  integra¬ 
tion  over  the  complex  shape.  However  this  is  also  not  the  way  of  using  simplex  integration. 

Giving  a  polygon  P1P2P2P4P3PG  with  P6  =  Pi,  such  that  Pt  rotate  at  the  same 
direction  as  from  ox  to  oy.  For  any  point  Pq,  the  algebraic  addition  of  the  2-d  simplex  vol¬ 
ume  (area)  P0P1P2,  P0P2P3,  P0P3P4,  P0P4P5  and  P0P5Pi  is  the  area  A  of  the  polygon. 
Let  P0  =  (0,0), 


1  5 


t=i 


1 

0 

0 

5 

1 

Xi  { 

X2  i 

1=1 

X\  i  X2  i 

1 

£l  i+1 

2-2  i+1 

Zl i+1  Z2 i+1 

In  Figure  2,  the  area  of  simplex  P0  P2  P3 .  Pq  P3  P4 ,  Pq  P4  P5  and  P0  P5  Pi  are  positive; 
the  area  of  simplex  P0P1P2  is  negative.  The  algebraic  sum  is  exactly  the  area  of  polygon 
P\  P2  P2  P 4  P5  P&  •  The  area  A  is  then  represented  by  the  coordinates  of  boundary  vertices. 

P5  P3 


Figure  2.  Addition  of  plus  or  minus  area  of  simplex 
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In  general,  simplex  integration  can  compute  ordinary  integrations  without  subdivid¬ 
ing  2-d  domains  to  triangles  and  3-d  volumes  to  tetrahedrons.  FEM  mesh  is  unnecessary 
for  simplex  integration.  Using  simplex  integration,  the  integration  of  any  n-dimensional 
polynomials  can  be  represented  by  the  coordinates  of  boundary  vertices  of  generally  shaped 
polyhedron. 


Two  Dimensional  Simplex  Integrations 

The  two  dimensional  integrations  here  in  this  section  are  used  for  two  dimensional 
DDA  algorithm.  Since  the  displacement  function  is  linear  function  of  coordinates  ( x,y ), 
the  integrands  are  of  degree  0,  1,2. 

A  2  dimensional  simplex  has  3  vertices 

Po,Pi,P2. 


Po : 

(®10, 

*2o) 

Pi  : 

(*11, 

*2l) 

P2  : 

(*12, 

*22 ) 

The  2  dimensional  coordinate  simplex  has  3  vertices 


U0,UUU2. 

U0:  (  0,  0) 

Ui:  (  1,  0) 

U2:  (  0,  1) 

The  following  coordinate  transformation 


(U1 1  U2)  (^i ,  x2) 


s  coordinate  simplex 


UolhUi 


to  normal  simplex 


P0P1P2 

(1-X)l"*)  +  *11  «1  +  x12  U2 

(1  —  53 1  «i)  +  X21  ul  x22  u2 
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x1 

Z2 


^10 

^20 


The  Jacobi  determinant  is 


j  _  ,  X2 ) 

D(ui,  u2) 

dx\  dx\ 

_  du i  du? 

du i  5u2 

_  *Ell  £10  #12  ^10 

^21  “  ^20  ^22  ~  ^20 

1  X10  X20 

=  1  x2i 
1  Z12  0:22 


Since  P0PiP2  is  a  2-dimensional  simplex  with  non-zero  volume,  J  is  non-zero. 
Translation  can  be  rewritten  as 


(uq, 111,112)  — >■  (1, £1,2:2) 


1 

= 

Uq 

+ 

Ui 

+ 

U2 

x\ 

=  X10 

u0 

+ 

xu 

Ul 

+ 

X\2 

U2 

X2 

=  X20 

Uo 

+ 

X21 

Ui 

+ 

X22 

U2 

Two  dimensional  simplex  integration 

Sp0plp2{rnl,m2) 


on  a  two  dimensional  simplex 


P0P1P2 


is  defined  as  normal  integration  times  the  sign  of  determinant  J. 


SpoP^Am,^) 

=  sign(J)  f  f  x™1x™2dx\dx2 

J  JP0PiP2 

=  sign(J)  /  f  x™lx™2 \ J\du\du2 

J  JUaUiUr, 


J 


Uo^U? 


'UqUiU? 

x™1  x™2du\du2 


Here 


2  m‘ 

*r  =  (£>/*«*) 

k=0 
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then  the  two  dimensional  simplex  integration  can  be  represented  by  the  following  basic 
forms  of  simplex  integration: 


Si  —  /  u^u^u^duidu^ 

J  Ju0UiV2 

y  yuo-f  +  «2==1 

—  I  I  UlQ  u'l  Ul2  dv,\du2 

J  j  Uo,U\}U2>0 

nl-~u2 

Uq°  du\du2 

_ 

rl  /  fl~U2 


Jo  {Jo 


u'l  (l  —  U2  —  Ui)  0 dill  dll2 


After  ii  times  of  integration  by  parts,  the  inner  integration  can  be  computed. 


/  11*1(1  —  112  —  ui)todui 

Jo 

yl— U2  1 

=  /  d(-^-r7tti1+1)(l-ti2-u,),° 

Jo  *1+1 

1  *1+1/1  \*0  |l~u2 

=  ^— 7U1  (1  —  t*2  —  «l)  0 

*1+1 

fl—U2  1 

-  ^—u\‘  +  1d((l-u2-uO'°) 

Jo  *1  +  1 

fl-u2  1 

=  -/  .  u!|1+1d((l  -  u2  ~  Hi)* 

Jo  *1+1 

=  /  T — — -Ui1+1(l  -  u2  -  **i)‘0  d 

Jo  *1+1 


-Ui1+1(l  —  u2  -  ui)'°  ldui 


.  t*l  +*0  —  1 


»  .  ,  1  1  V-  *  */  A 

7o  ^i  +  l 

: J  («,  +(i°)(i11+  2j“"+2(1 "  “2  “  Ul)‘°  2<i“1 

=  ■ofa-1)(i°-2)  u‘,+3(1  -  u2  -  ui)'°-3dui 

Jo  (*i  +  l)(*i  +  2)(*i  +3) 

/1_U2  *o(*o  ~  l)(*o  —  2) ...  2  ^1+io_ 

7o  (*1  +  l)(*i  +  2)(*i  +  3) . . .  (*1  +  *o  —  1)  1 

(1  —  112  —  u\)1  du\ 

:  C~U2  *o(*Q  —  1)(*0  —  2)  ...  1  ii+io 

Jo  (*1  +  l)(*i  +  2)(*i  +  3) . . .  (*i  +  *0)  1 

=  f~"  7j£r!v-i1+,0^i 

Jo  (*1  +  *0)! 
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1 


1  —  u2 


Z0!2l! 


(ii  +  io  +  1)! 


d(u\i+l0+1) 


(io  +  it  +  1)! 


io+ii  +  l  |1_U2 


i0!ii! 


( i\  +  *o  +  1)! 


(1  ~“2) 


H+io+1 


rl  f 1-U2  .  j 

5!=  I  u%du2  /  <(1  -u2  -wi)  °dui 


=  f1  u^(l-u2)io+il+1du, 

Jo 


i0lix\ 


(io  +  it  +  !)• 


ioHi!^1- 


(io  +  *i  +  *2  +  2)! 


Based  on  the  previous  formula  of  Si ,  the  two  dimensional  integration  of  polynomi¬ 
als  of  degree  0,  1,  2  can  be  computed. 


SpoPlP2(0,0)=sign(J)  f  f  dXldx 

J  J  Pn  Pi  P2 


=  J 


Pq  P\  P'2 
du\du2 


UoUx  U 2 


5PoPlP2(0,0)  =  sign(J)  /  f  dxxdx2 

J  JPoPlP? 


=  J- 


0! 


P0P1P2 


(0  +  2)!  2 

5poPip2(1,0)  ^stgn(J)  /  I  x1dx1dx2 

J  J  Pq  Pi  P2 

=  J  /  (xio^o  +  *11^1  +  X\2U2)du\dU'} 

J  JU0U1U2 

f  (xiO  +  *11  +  *12) 
x\dx\dx2 


=  J- 


1! 


(1+2)! 

5pop1p2(1,0)  =  sign(J)  J  J ^ 


P0P1P2 


1 


=  7^(*io  +  *11  +  *12) 
6 


5p0P1P2(°>  !)  =  sign(J)  f  j  x2dxxdx2 

J  JP0P1P2 
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—  20  X21  +  ^22) 


^PoPi P2 (2,  0)  =  sign(J)  /  /  x\dxxdx2 

J  JPqPxPi 


=JIL 


UoUiUi 


(xiouo  +  inui  +  X12U2)2  duidu 


=  J 
+  J 


2! 


(2  +  2)! 

1! 


(*10  +  +  *^12 ) 


(2  +  2)! 

(0  +xioxn  +xi0a;i2 
+2:11X10  +0  +X11X12 

+X12X10  +X12X11  +0) 


Spop1pi(2,0)  =  sign(J)  /  /  x\dx\dx2 

J  J  Pq  Pi  P^ 


4 


(2xioXio  +X10X11  +X10X12 

+X11X10  +2xijXh  +X11X12 
+x12x10  +X12X11  +2x12X12) 


SpoPj, p2(0, 2)  =  sign(J)  /  /  x\dx\dx2 

J  J  PqP  1  P2 


=  J 
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(2x2ox2o  +X20X21  +X20X22 

+X21X20  +2x2ix2i  +X21X22 
+X22X20  +X22X21  +2x22x22) 


Sp0p1p2(l,l)  =  sign(J)  /  /  xix2dx\dx2 

J  JP0P1P2 

=jJL 


(xioUo  +  X11U1  +  X12U2) 

'UoUiUu 

(x2qUq  +  X21U1  +  x22u2)duidu2 
2! 

=  J—  “(X10X20  +X11X21  +X12X22) 


+  J 


(2  +  2)! 

1! 


(2  +  2)! 

(0  +^10^21  +#10^22 

+^11^20  +0  +X11X22 

+  ^12^20  4”  X 1 2  X  2 1  +0) 


5'p0PiP2(l,  1)  =  sign(J)  /  /  x\x2dx\dx 

J  JP0P1P2 


=  —  J(2xio*20  +  2X11*21  +  2X12X22 

+  *10*21  +  *10*22  +  *11*20  +  *11*22  +  *12*20  +  *12*2l) 


(2xio*20 

+*11*20 

+*12*20 


+*10*21  +*10*22 
+  2xnx2l  +*11*22 
+*12*21  +2*12*22) 


Three  Dimensional  Simplex  Integrations 

The  three  dimensional  integrations  here  in  this  section  are  used  for  three  dimensional 
DDA  algorithm.  Since  the  displacement  function  is  linear  function  of  coordinates  (x,  y,  z ), 
the  integrands  are  of  degree  0, 1,2. 

A  3  dimensional  simplex  has  4  vertices 

P0,Pl,P2,Pz. 

P 0  •  (*10,  *20)  *30 ) 

Pi  ■  (*11,  *21,  *31 ) 

P2  :  (*12,  *22,  *32) 

P3  ■  (*13,  *23,  *33) 

The  3  dimensional  coordinate  simplex  has  4  vertices 

U0,UUU2,U3. 


Uq  : 

( 

0, 

0, 

0) 

Ui  : 

( 

1, 

0, 

0) 

U2  : 

( 

0, 

1, 

0) 

U3  : 

( 

0, 

0, 

1) 

The  following  coordinate  transformation 


(111,112,113)  -»  (*1 , *2 , *3 ) 


s  coordinate  simplex 


UoUiUiUi 


to  normal  simplex 


P0P1P2  Pi¬ 


rn 


X1 

=  Xio 

(i 

—  Ei  u»)+ 

xn 

U1 

+ 

X12 

U2 

+ 

^13^3 

X2 

—  X20 

(i 

-  Ei  «.•)+ 

X21 

U\ 

+ 

X22 

U2 

+ 

^23^3 

X3 

=  ^30 

(i 

-  Ei  u*)+ 

X31 

U\ 

+ 

X32 

u2 

+ 

X33U3 

The  Jacobi  determinant  is 


J  _  Djxu  %2,  X3) 

D(u  1,  u2,  uz) 


dx\ 

dx\ 

dx  1 

dui 

du2 

du  3 

dx2 

dx2 

dx2 

du\ 

du2 

duz 

dxs 

dxz 

dx  3 

i 

du  1 

du2 

du3 

xn 

~  Xi0 

Xl2  ~ 

-  X10 

^13 

-  XlO 

= 

X21 

-  ^20 

X22  - 

-  X20 

^23 

-  X20 

X31 

“  £30 

X32  - 

-  2:30 

Z33 

~  Z30 

1  *^10  ^20  Xzo 

_  1  Xll  X21  Xzi 

1  Xi2  X22  XZ2 
1  2:13  X23  Xzz 

Since  PqPiP2Pz  is  a  3-dimensional  simplex  with  non-zero  volume,  J  is  non-zero. 
Translation  can  be  rewritten  as 


(«0, 

«1, 

) 

(1, 

xnx2, 

,X3) 

1 

= 

Uq 

+ 

U\ 

+ 

u2 

+ 

U3 

Xi 

=  x10 

Uq 

+ 

xn 

Ul 

+ 

X\2 

U2 

+ 

Xl3 

U3 

X2 

=  X20 

n0 

+ 

X21 

Ul 

+ 

X22 

U2 

+ 

X23 

U3 

X3 

=  X30 

no 

+ 

X31 

U\ 

+ 

X32 

u2 

+ 

X33 

U3 

Three  dimensional  simplex  integration 


Sp0plp2p3{mum2,mz) 


on  a  n  dimensional  simplex 


PqP\P>P:\ 


is  defined  as  normal  integration  times  the  sign  of  determinant  J. 


Sp0PiP2P3(mUm2,m3) 


=  sign(J)  J  J  Jp 
=  sign(J)  J  J  J 


P1P2P3 

UqUi  U2U3 


x™1  x™2x™3  dx\dx2dxz 
x™1  xpx?*  \J\duidu2duz 


=  J 


'u0ul  U2U3 


x™1  x™2  x™3  duidii2duz 
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Here 


*r = cy:  xikUk) 

k= 0 

then  the  three  dimensional  simplex  integration  can  be  represented  by  the  following  basic 
forms  of  simplex  integration: 


*-/// 

J  J  JUoUXU2 

r  r  /k«o+«i 

J  J  J  Uo,Ui,U' 

nl  —  u3  r 

Jo 


/  /  ul?u\1ut^ut^duidu2du3 

J  JUoU^Us  ‘  ' 

/  /  ulQUl^ul2Ut^du\du2duz 

J  Juo,U\ ,U2,Uz'>0 

1  rl—Uz  rl~U3“w2 

I  I  u^u'^Ui  ^o° duidu2(iu3 


/*1  /*1  —  U3  /  rl—us  —  u^  .  \ 

=  J  J  u£ul2  yj  u\  (1  —  u3  -  U2  -  Ui)l°dui  J  du2du3 

After  ?i  times  of  integration  by  parts,  the  inner  integration  can  be  computed. 


*1—  Us  —  U2 

u\l  (1  —  uz  —  U2  —  ui)l°du\ 

) 

f*l— U3-U2  1 

d(- - r«\1+1)(l  -U3-U2-  Ul)t0 

)  *1+1 

1  ,11  +  1/1  \*0  1 1  «3  1*2 


Wj  (1  “  Uz  -  U2  -  **l) 


*1  + 1  1  v 
/'1-US-U2  1 


*1+1 


tjl  +  1C?((  1  -U3-U2-  «l)‘°) 


-i 

■r 


1—  U3  —  U2 


*1  +  1 


tj1  +  1d((l-U3-U2-Wl)l°) 


1— W3  — 


-wl11+3(l  —  u3  —  u2  —  ui)t0  dui 


/  •  ?°  T  ^11+1(1  -U3-U2-  Ui)t0  ldu\ 

Jo  *1  +  1 

/•1-“3-U2  *o(*o-l)  il+2/1  ^0-2^ 

L  (M  +  !)(,.+  2)-  (1-"3-“2-”l)  dUl 

=  ■°c°-1)(r2)  x+3a «,)' 

Jo  (*l  +  1 )  ( *  1  +  2)(*i  +  3) 

;  ~‘a  »o(-o-l)(<o-2)...2  ui.+i.-i 

Jo  (*1  +  l)(*i  +  2)(*1  +  3) . . .  (ii  +  io  —  1) 

(1  —  1/3  -  U2  —  U\f  du\ 

,  rua~u2  *o(*o  ~  i)(*o — 2) ...  1  n,1+i0(ju 

Jo  (*i  +  l)(*i  +  2)(*i  +  3) . . .  (ii  +  *o)  1 


*  1  — *  U3  —  U2 


-1—  W3  —  U2 


r*l—  W3  —  U2 


i’l  +  io-l 
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*1  u3  u%  *  I*  I 

)  (*1  +  *o)‘ 


*1~U3-“2  *0!*i ! 


d(«j1+<0+1) 


(it  +  *o  +  1)! 


io+!l  +  l  |l-u3-«2 

:(io  +  ii  +  l)!'  1  io 


(*i  +  *o  +  1)! 


(1  -  U3  -  U2) 


*l+*0+l 


1  /*!— 1*3 


•1— u3  — U2 


S[  =  I  I  u'zulfduiduz  I  u\: 

Jo  Jo  Jo 

r  1  pl—Uz 

-  /  *3  /  „>2n  _  _  „.V0+'»+1 


^(1  -U3-U2-  ui)todiii 


«22(1  -  «3  -  ^2) 


du2du3 


lo'-ll- 


1 0  Jo 
/•l 


(i0  +  *i  +  1)! 


=  /  u*3(l  -  u3),0+il+,'2+2rfu3p - 

7o  (*o  +  *i  +22  +  2)! 


i0!ii!i2!i3! 


*3  +  3)! 


Based  on  the  previous  formula  of  Si,  the  three  dimensional  integration  of  polyno¬ 
mials  of  degree  0,  1,  2  can  be  computed. 


‘S,fbPlftPa(M»0)  =  sign(J) 


dx 1 dx2dx 3 


PO  Pi  P’2.  P 3 


(0  +  3)!  6 


^PoPiPaPsO^O)  =  sign(J) 


x\dx\dx2dx3 


Pq  Pi  P’2  P3 


(a^io^o  +  a?iiwi  +  ^12^2  +  £13^3) 


J  J  JUqUxUvUz 
du\du2du3 

=  J ^  +  Xl1  +  Xl2  °°13 ) 

SpoPi P2P3 (19  0, 0)  =  sign(J)  /  /  /  xidxxdx2dx3 

J  J  J PoPiP2^3 

—  «/—  (#io  +  ^11  +  ^12  +  ^13) 
5,PoPiP2/:>3(0j  1?  0)  =  62 <771  ( J)  /  /  /  x2dx\dx2dx3 

J  J  JPoPiP^Pz 
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=  J —  (£20  +  x 21  +  x 22  +  x 23) 


‘S,poPiP2p3(0, 0, 1)  —  sigvi^J')  III  x3dx\dx2dx2 

J  J  JPqPi  P2P3 


J ckA  (2-30  +  2-31  +  £32  ”1“  £33) 
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S'poPiP2p3(2,0,0)  =  5^(7)  /  /  /  x\dx3dx2dx3 

J  J  ./P0P1P2P3 


J 


in 


(x10UO  +  XllUl  +  X12U2  +  Xi3U3y 


UoUi  U2U3 


du\du2duz 
=/ 2! 

+  J 


(2  +  3)! 
1! 


(xio  +  xn  +  *12  +  ^n) 


(2  +  3)! 

(0  +£io£n  +£io£i2  +£io£i3 

+£ll£lO  +0  +£ll£l2  +£ll£l3 


+£12X10 

+Zi2£ii 

+0 

4*2:12X13 

+^13^10 

+2:132:11 

+2:132:12 

+0) 

5'pop1p2p3(2,0,0)  =  sign(J) 

Ijl  x\dx\dx2dx3 

J  J  JP0P1P2P3 

=  J— 
120 

(2xio^io 

+XioXn 

+2:102:12 

+£102:13 

+X11X10 

+2xiixn 

+2:11X12 

+£n£i3 

-\-x12X10 

+2:122:11 

+  2X12X12 

+£12£13 

+  ^13^10 

+2:132:11 

+X13X12 

+2X13X13) 

5'poP1P2P3(0,2,0)  =  sign(J) 

/  /  /  x\dx\dx2dx$ 

J  J  JP0P1P2P3 

=  J— 
120 

(2x202-20 

+£202:21 

+2:202:22 

+X2o£23 

+X21X20 

+2x2ix2i 

+2:212:22 

+  £2l£23 

+Z22Z20 

+2:222:21 

+2x222:22 

+X22£23 

+£23  £20 

+2:232:21 

+2:232:22 

+ 2x23£23 ) 

‘S,poPiP2P3(0’0>2)  =  sign(J) 

I  x\dx\dx2dx$ 

J  J  j  Po  Pi  P2  P3 

=  J 


120 
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(2X30^30  +^30^31  +^30^32  +^30^33 

+^31^30  +2X31X31  +X31X32  +  0:31X33 

+  ^32^30  +^32^31  +  2x322-32  ~hX32X33 

+^33^30  +^33^31  +^33^32  +2x33^33) 


SpoPlp2p3(l,  1,0)  =  sign(J)  /  /  /  xix2dx1dx2dx3 

J  J  JP0P1P2P3 

=  d  I  I  I  (xiotio  +  X11U1  +  X12U2  +  X13U3) 

J  J  Ju0UiU2U3 

(X20U0  +  X21U1  +  X22U2  +  X23U3)dU\du2du3 

2! 

=  J-( 2  ^'^(^10^20  +  XUX21  +  X12X22  +  0:13X23) 

11 

+  J(2  +  3)! 

(0  +0:10X21  +0:100:22  +0^100:23 


+X11X20  +0 

+X12X20  +0:120:21 

+X13X20  +0:13X21 


+  X11X22  +X  110:23 
+0  +0:120:23 

+0:130:22  +0) 


Sp0p1p2p3(1,  1,0)  =  sign(J)  /  /  /  x1x2dx1dx2dx3 

J  J  JP0P1P2P3 


J 


120 


Sp0p1p7p3( 


(2^10^20 

+0:100:21 

+0:100:22 

+  Xio^23 

+  X11X20 

+2xiix2i 

+0:11X22 

+  ^11^23 

+Xi2^20 

+X12X21 

+2X12X22 

+X12X23 

+  X13X2O 

+0:130:21 

+  X13X22 

+  2X13X23) 

=  sign(J) 

l  i  jp0Pi 

x\x3dx\dx2dx3 

P2P3 

=  j— 

120 

(2xi0x30 

+X10X31 

+X10X32 

+^10^33 

+  XHX30 

+2xnX3i 

+0:11x32 

+^11^33 

+^12^30 

+0:120:31 

+2xi2x32 

+  Xi2^33 

+  ^13^30 

+0:130:31 

+X13X32 

+2X13X33  ) 

=  sign(J) 

/  /  /  x2x3dxidx2dx3 

J  J  JPoPiPiPs 

=  j— 

120 

(2^20  ^30 

+X2oX3i 

+  X20^32 

+ x2qx33 

+X21  X30 

+2X21X31 

+^21^32 

+0:210:33 

+  ^22^30 

+  X22  0:31 

+  2X22^32 

+  X22X33 

+  ^23^30 

+0:230:31 

+  ^23^32 

+  2X230:33) 
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Three  Dimensional  Simplex  Integrations 
to  the  Coordinate  Simplex 

A  n  dimensional  simplex  has  n+1  vertices 


Po,Pi,P2,  •  •  •  ,Pn- 


Po : 

(*10, 

*20, 

2-30 ,  •  • 

•  ?  *^rxO  ) 

Pi  : 

(*n» 

*21, 

*31,  •• 

•  ?  ^nl  ) 

P2: 

(*12, 

*22, 

*32,  •• 

•  5  ^n2  ) 

Pn  : 

(•Tin) 

*2  n  > 

*3 n ,  *  • 

*  5  *^nn) 

The  n  dimensional  coordinate  simplex  has  n+1  vertices 


0b, 

uU2, 

,  •  • • ,  Un  • 

U0: 

( 

0, 

0, 

0,  ..., 

0 

) 

Ui  : 

( 

1, 

0, 

0,  ..., 

0 

) 

U2  : 

( 

0, 

1, 

0, 

0 

) 

( 

0, 

0, 

0, 

1 

) 

The  following  coordinate  transformation 


(1) 


(til ,  W2 ,  ^3 ,  •  •  •  ,  ^n)  ^  (xi,3/2,*35*  •  •  ?  *n  ) 


transfers  coordinate  simplex 


U0UiU2...Un 


to  normal  simplex 


*1 

=  *10 

(1- 

*2 

=  *20 

(1- 

*3 

=  *30 

(1- 

*n 

=  *n0 

(1- 

The  Jacobi  determinant  is 


P0P1P2.. 

■Pn 

e: 

«,)+ 

*11 

Ux 

+ 

Ei 

wt)+ 

*21 

Ui 

+ 

n 

2^1 

Wi)+ 

*31 

Ui 

+ 

V-^\7l 

2^i 

l*t)+ 

*nl 

Ui 

+ 

*12 

U  2  +  .  . 

.  + 

%1  n 

Un 

*22 

U2  +  .  . 

.  + 

n 

*32 

U2  +  .. 

.  + 

^3  n 

*n2 

u2  +  ... 

.  + 

^nn 

(2) 


j  D{x\,  X2,  *3,  X7 j) 

D(uu  u2,  u3,  ...,  «„) 
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dx\ 

dx\ 

dx  1 

dx  1 

du\ 

du2 

du3 

<9Un 

dx2 

dx2 

dx  2 

dx2 

du  1 

dU2 

du  3 

dun 

dx3 

dx3 

<9  3:3 

dx  3 

du\ 

du2 

du  3 

5un 

dxn 

dxn 

<9a:n 

C?£n 

du  1 

du  2 

<9«3 

5un 

xn  - 

XlO 

£12  “ 

£10 

£13 

—  X10 

•  •  £ln 

~  £10 

£21  ~ 

X20 

X22  - 

£20 

£23 

—  X20 

£2  n 

”  £20 

= 

£31  — 

X30 

X32  ~ 

£30 

£33 

-  X30  • 

•  •  £3n 

~  £30 

£ nl 

XnO 

Xn2 

£n0 

£  n3 

XnO 

•  ■  xnn 

XnO 

xn  - 

XlO 

X21  - 

£20 

£31  ■ 

-  X30  • • 

•  Xnl 

—  xno 

xu  ~ 

XlO 

X22  - 

£20 

£32  • 

-  X30  .  . 

■  Xn2 

—  xno 

= 

xn  ~ 

XlO 

^23  “ 

£20 

£33 

-  X30  • • 

•  £n3 

£n0 

X\ n 

XlO 

£2  n 

£20 

£3n 

-  X30  •  • 

•  Xnn 

”  £n0 

1 

x\o 

£20 

£30  •  « 

-  •  XjiO 

1 

«n 

X21 

£31  •  < 

■  .  £nl 

1 

X12 

X22 

£32  '  ' 

.  .  £n2 

1 

Xin 

X2  n 

£3n 

>  ♦  xnn 

Since  P0-P1-P2  . . .  Pn  is  a  n-dimensional  simplex  with  non-zero  volume,  J  is  non-zero. 
Translation  (2)  can  be  rewritten  as 


(«0 

,«2, 

u3?... 

yUn 

)-H 

(Ml, 

,£2, 

X3  ,  •  •  • 

•  1  Xn) 

1 

= 

Uo 

+ 

Ul 

+ 

U2 

■  + 

un 

£l 

=  £10 

Uo 

+ 

£11 

Ul 

+ 

£12 

U2 

.  + 

£ln 

Un 

£2 

=  £20 

Uo 

+ 

£21 

Ul 

+ 

£22 

U2 

.  + 

£2n 

Un 

£3 

=  £30 

Uo 

+ 

£31 

Ul 

+ 

£32 

U2 

.  + 

£3n 

Un 

Xn 

“  XnQ 

Uo 

+ 

£nl 

U\ 

+ 

£n2 

U2 

+  •  • 

•  + 

Xnn 

Un 

/  1  \ 

( 1 

1 

1 

1  ^ 

(U0\ 

£l 

£10 

Xn 

£12  ■  - 

■  •  £ln 

U 1 

£2 

£20 

2-21 

£22  •  • 

•  •  £2n 

«2 

£3 

£30 

X$1 

£32  •  - 

•  •  £3n 

u3 

\xn  / 

Vx„o 

£nl 

£n2  •  ■ 

Xnn  / 

\Un  / 
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The  determinant  is 


1 

1 

1 

1 

1 

a-'io 

%20 

X30 

XnQ 

Z10 

xn 

Xi2  . 

•  •  ^1  n 

1 

xn 

X21 

X31  ■ 

•  *  %nl 

X20 

X21 

X22  ■ 

•  *  ^2  n 

1 

Z12 

X22 

X32  ■ 

•  *  *^n2 

X3Q 

X31 

X32 

•  •  ^3n 

1 

®13 

X23 

X33  ■ 

•  *  ^n3 

^nO 

%nl 

xn2 

•  *  ^nn 

1 

X\n 

3'2n 

X3  n 

■  •  ^nn 

From  formula  (4)  (u o,  •  •  • ,  «n)  can  be  computed 


/U0\ 

( 1 

1 

1  . 

•'  1  ^ 

-1 

/  1  \ 

Ul 

^10 

xn 

xn  . 

.  .  X\ n 

Xi 

U2 

X20 

Z21 

x22  • 

.  .  %2n 

x2 

Us 

X30 

2:31 

^32  • 

■  •  ^3  n 

X3 

\un/ 

V  ^n0 

%n  1 

Xn2 

•  •  ^nn  / 

V  Xji  J 

Definition  of  N  Dimensional  Simplex  Integration 


Simplex  integration 


Sp0p1p2...pn{r)-i1,m2,m3, . . .  ,mn) 

on  a  n  dimensional  simplex 

P0P1P2  •  ■  ■  Pn 

is  defined  as  normal  integration  times  the  sign  of  determinant  J. 

Sp0p1p2...pn{mi,m2,mz,. . .  ,mn) 


=  sign(J)  f  f  I  ■  f 

J  J  J  JPoPl 

=  sign(J)  Iff 


P2-Pn 


x™1  x™2x™3  ■  •  •  x™n  dx\dx2dxz  ■  ■ .  dx? 


x™1  x™2x™3  ■  ■  •  x™n  \  J\duidu2du3  . . .  du , 


'UoU^-Un 


J 


x™1  x™2 x™3  •  •  ■  x™n  du\du2duz  . . .  dun 


Ju0ulu2...un 

In  order  to  compute  the  integration,  coefficients  of  invariants  u\  in 


^mi 


71 

(^2 xikUk ) 


of  formula  (6),  the  following  formula  (7)  is  useful. 


The  formula 


<X»  =£ 


io+ii+is+  - +in=n»  (j0  +  ji  +  J2  +  •  •  •  +  in)! 


joiil  list  •••(in  >0 


jo!ji!i2!  •  •  "in! 


'0°alla22  •••<"  (7) 


can  be  found  in  algebra  books,  however  a  brief  derivation  are  given  here  for  convenience. 
In  case  of  n=0,  formula  (7)  is  correct: 


<£“>  =£:::!“»=< 
k= 0  “  J 


In  case  of  n=l,  formula  (7)  is  correct: 


/V-  \  _  v-io+J1  =  m  (io  +  Jl)! -io-ii 

^  ^  ^  io.ii  >o  io  !ii !  0  1 

Assuming  formula  (7)  is  correct  for  n: 

/V'  x  y-v7o+ii+is  +  ---+in  =  m  (j0  +  jj  +  j2  +  •  •  •  +  in)!  ,  ,  ,2  n 

/c=0 


then  the  following  computation  shows  formula  (7)  is  correct  for  n+1 


(£«*) 


=  ((^  Ctk)  +  On+l) 


Ein+jr.+1-m  (zn  +in+l)!  ,r^  \  in  +  l 

*n,in  +  l>0  t„!i„+l!  "+1 

_  ^-^*n+jn  +  l  =  ni  (in  +  jn+1)! 
in  ,jn  +1^0 

v^io+ii+j2+.»+i„=*n  (jo  +  ji  +  h  +  •  •  ■  +  in)!  ^  Jn+1 

/  .  .  .  .  .  .  _  •  I  '  I  "  I  FF  a°  G1  a2  *  *  *  an  an+l 

Ein  “1“  jn  4-1  =  ni  1 \JO  Jl  H“  J2  "4"  •  ■  jn 

in  ,jn  +  l  >0  '  jo  Jl  j  J2  i  *  •  •  jjn  >0 

(^n  +  jn+l)-1  (jo  +  jl  +  J2  +  •  •  •  +  Jn)!  ;0  j,  j2  _  Jn  Jn  +  1 

* - r-p - j - - FI  a0  al  a2  "  *  an  an+l 

^n.Jn+l*  j0*jl-j2 . Jn* 
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V*n+in  +  l  =  m  V^JO+il+i2  +  ...+in=*n 


'*„,Jn  +  l>0  " 

(in  H“  Jn+l)! 


Oo.il  ,i2v*.in>0 


V  n  1  Jo  ?!  J2 

io!ii!i2!...Win+i!  °  1  2 

jo  +il  4- J2  +  . .  .+in  +in  + 1  =  m 


,  0in  a-?n  +  1 
an  an+ 1 


*-^jo,jl  ,J2  ."-On  ,in  +  l  >° 

(jo  +  jl  +  J2  +  «  ■  •  +  in  +  Jn-f-1 ) » 


<*0°  al1  a2 


„jn  „•?"  +  ! 
°n  an+l 


Denote 


from  the  previous  algebraic  formula, 


-i-E' 


E*tO  i  1  " 
lin.tn  A 


■T  =  £>*,*) 


ilo  +  ili+il2+—+iln=rni  (i/o  +  i/l  +  il2  +  •  •  •  +  i/n)! 


*10  *11  *12  .  .  .  7/*ln  *10  *11  *12  ...  r*ln 

a0  Uj  u2  x/0  x/2  xln 


Substituting  (8)  into  (6),  integral  (9)  is  obtained. 


Sp0p1p2...pn(mU7n2,m3,...,mn)  =  J  /  /■••/ 

J  J  J  JUoUiU2...Un 

^*104**11 +*12+ •.•+*!« =*^1  (z’io  +  in  +  ii2  +  .  • .  +  im)! 

*10)*11  J*12»--M*ln  >0  lin!lii!tio! 


/  v*204-*2l4"*224-. 

/  *20, *21  .*22 


iio!tnUi2!  •  •  -!*in! 

«#*107,*11  7 ,  *  1 2  .  .  .  7|*ln  ™*10  ^*11  12  .  .  .  r*ln  ' 

u0  al  a2  an  x10  X11  x12  X1  n  , 


.+i2n  =  m2  (l20  +  l2l  +  *22  +  •  •  •  +  *2n)! 


*20**21**22*  •  •  -l2n! 


*207/*217/*22  .  .  .  7/*2n  ^*20^*21  ^,*22  .  .  .  „*2n  \ 

a0  U1  U2  un  x20  x21  x22  ^2n  ) 


^^*3o4-*3l4-*324-...4-*3n  =  *n3  (i30  +  i31  +  i32  +  .  .  .  +  23n)! 

-/*3O,*3i,*32,-.-,*3n>0  *30^31^32!  •  •  "i3n! 


7/*307/*317/*32  .  .  .  7  .  *  3  n  ~*30  _*31  *32  .  .  .  „*3n  \ 

a0  al  u2  an  x30  ^31  x32  x3n  ) 


,«'n0+»nl+«n2+...+»nn  =  mn  (inQ  +  jnl  +  in2  +  .  .  .  +  tnn)! 


^  *  n  0 .  *  n  1  j  *  n  2 1  •  •  • » *  n  n  ^  0 


du\du2du^  . . .  dun 


^nO*^nl'in2*  *  •  ••^nn* 

7/*n07,*nl7/*n2  .  .  .  9/*nn  TlnOT *  n  1  ~  *  n  2  .  .  .  7,lnn  \ 

u0  al  u2  un  xnO  xn  1  xn2  xnn  J 
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where 


Sp0p1p2...pn(m1,m2,m3,..  .,mn) 


='E 


*10  +  *11  +  *I2+ •••+*!  n”  ”"11 


\*30”H31  +  *32+  ••  •  +  *3n  =  *™3 
'*30  ,*31  i*32>“-i*3n  >0 


E*20  +  *21"H22+...-H2n  — *™2 

ion.ioi  .in . >0 


-  E 


*20  i*21  )  *22  j  •  •  *  >  *2n  >0 
2n0  +  inl+ln2+-"  +  *nn=:r^r 


*rtO  ?*nl  >*n2j***)*nn 


m\\ 

m2! 

*io!*n  1*12!  •  •  - 

Un- 

*20-*21  ^22?  • 

*  *^2n* 

m3! 

m„! 

*30^31  ^32!  •  •  •! 

&  3  n J 

^nO  *^nl  *^n2  •  •  • 

•^nn  ■ 

TJior*n  ~H2  . 

x10  X11  x12 

*  1  n 

xln 

*20  ™*2i  *22 

x20  '*'21  x22 

.  .  .  ^r*2n 
x2n 

*30  rp  *31  /y.  *32  .  , 

X30  ^31  ^32 

.  'Tt3n 

x3n 

T*n0  r*nl  r^n2  . 

xn0  xnl  xn2 

.  ,  ry»*nn 

^nn 

fff-f 

«0  ullu2a  ■  •  •  un 

du\du2 

WQVlUi...Un 


il  —  ill  +  *2/  +  *3/  +  ‘  ‘  +  *n(i  ^  —  0,1,  2, 


Formula  of  N  Dimensional  Simplex  Integration 

The  coordinate  simplex  integration  in  (10)  is  another  form  of  well  known  Dirichlet 
integration  in  classical  analysis.  The  direct  computation  of  the  integration  is  given  below. 


-///-/„ 

=///•/; 

nl-«n  /O 

Jo 

ln  *n  —  1  *n~2 


/  /  *••  /  *  •  *  U^ncfoi(iu2<iu3  •  •  *  dun 

J  J  JUoU1U*...Un 

r  r  /'Uo~f-ui4'u2'f..-+Wn:::::l 

/  /  •••  /  Uq’u^Uj2  ■  ■  •  ul^du\du2du3  . . .  du, 

J  J  J«0)Ui,«2r-tun>0 

1  /*  1  ti  yi  /*  1  ri  n  1  /*  1  V>n  **  n  —  1  **2 


U^n  un"_i  «„_2  ’ ' '  u\1ut0°duidu2du3  . . .  duT 


nl— U„  rl—  un  —  «n_i  /-1-Un-tin-l - « 3 

Jo  JO 


/  /»!  —  un—  un_i - u  2 


l1  ~  1  Ufc)  )  ^2^3  •  *  •  du 


After  ii  times  of  integration  by  parts,  the  inner  integration  of  (1 1)  can  be  computed. 


fl-Un-Un_i - U  2 


/  ^  n  \  *  o 

H1~J2k=iUk)  dui 
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r 1  u"  "“-1  "■  “2  1  .  ,  ,  /  y — ,n  \io 

=J0  )(1~Et=,u*) 

=i7TT'‘i,+,(1-El=1“‘)iolo"“""‘ 

/•l-un-un_! - u2  j  /  n  \«o 

~l  iT TI"«  d((i-Efel“‘)  ) 

- irh“'',+1  (» -  EL,  «*) 

=  C - 

,/o  (®i  +  l)(*i  +  2)  V  ^—sk-i  ) 

_  f1— *0(*0  -  1)(*0  -  2)  ii+aA  V0-3j 

~J0  (>l  +  l)(ii  +2)(«!  +3)  1  (‘  Etal“*)  d“' 

=  f1—-— »  i0(,0-l)(,o-2)...2 

Jo  (*i  +  l)(^i  +  2)(ii  +  3) . . .  (ii  +  zo  —  1)  1 

(* -E«  •*)■*- 

=  /*— - -  M<o-D(.o-2)...l  <1+io, 

7o  (*i  +  l)(*i  +  2)(t’i  +  3) . . .  (ii  +  ?o)  1 

-  f1  *"  Un_1  “2_*oM_  ,1+<0, 

7o  (*i+*o)!  1 

=  /  <0,>1< _ d(V1+,0+1') 

7o  (*i  +  *o  +  1)-  1  j 

*0!*l!  io-f-ij  -f-l  ll-Un-Un-X - U2 

“(i0+ii+l)!Ul  lo 

«o!*i !  (  'r-^n  \»i+*o+i 

~(*1  +*0  +  1)!  V1  ~  ^ k=2  Uk) 

/•l  rl-Un  /*l-Un-«n-l  /*1-Un-U„_i - «3 

Si=  /  /  •••/ 

JO  JO  JO  JO 

- “2  .  ,  y0 

j  u\1\^~/^k  ukj  duidu2<iuz  . . .  dun 

nl-Un  rl—Un -Un-1  /*1  W  n  —  1 - «4 

Jo  Jo 

r\-un~un.l  «3  .  „  Jo+il  +  1 

du^du^dui  . . .  dun— - EE - — 

(*o  +  *i  +  1)! 


(12) 
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1  rl  —Un  rl  —  Un—  Un_i  pi— Un  U  n  —  1  ***  U  4 


nl  — un  /* 

Jo 


r *o  +  *i  +  *2+2 

(1  —  2_^  uk)  du^du^du^  . .  .  du 


7/in,*n-l  *n-2  .  ?/*4  *3 

un-l  Un- 2  U4  u3 


t’o!ti!t2! 


~  Z^fc=3  Uk)  "  (*o  +  *!  +  t2  +  2)! 

1  r  1  —  U  n  r  1  Un  tXn  —  1  /*l-1in-Un_i  «5 


rl  rl  —  «n  ri 

7o  Jo  Jo 


7/*n?/tn“1?iln“2  •  .  .  7/U 

un  Un-1  Wn-2  U4 


rl  «n  u  n  —  1  "*  **4  — *0  +  *l  +  *2  +  2 

JO 

,  ,  ,  ,  *o!*i!*2! 

auiduidus  . . .  dun—  ■  :  ■  :  p7vT7 

(to  +  *1  +  *2  +  2)! 

rl  rl  un  *  -  ,  •  ,  ,  • 


n;  \*0  +  *l  +  *2+‘*’“Hn-2  + ™—2  r  j 

unun-l  (1  —  un-l  —  un)  <tUn-.\aUr 


i0!ii!i2!  •  •  .!in-a! 


(to  +  *1  +  *2  +  •  •  •  +  t'n-2  +  n  2)! 

:  f 1  uj,"(l  -  «n)*0+<1+,2+  +*"_1+n_1^u» 

Jo 


to!ii!i2!  •  •  -!tn-i! 


(to  +  *i  +  *2  +  •  •  •  +  tn-i  +  n  —  1)! 

_ _ t'olz'i  !i2!  •  •  -!in! _ 

(to  +  *1  +  *2  4"  •  •  •  +  in  +  ft)! 

Substituting  equation  (13)  into  (11),  the  integrals  of  a  n  dimensional  simplex 

Sp0plp2...pn (mi ,  m2,  m3, . . . ,  mn) 

is  obtained. 


Sp0p1p2...p„  (mi,  m2,  m3, . . .  ,mn) 


JE 


*lo4_*ll'4**12~K**“Hln  —  ml 


'*10  >*11  ,*12i-->*l 

^*30  +  *31  "f  *32"K”H“*3n:=*T*3 

'*30  >*31  > *32  >"*>*3n  ^0 


E*20  +  *21+*22  +  .--+*2n  — ™2 

l*>n  .loi  .Ion . ton  >0 


-  E 


*20  ,*21  ,*22  ,  •  •  *,*2n  >0 
*nO  +  *nl  +  *n2  +  -+*nn  =  ”*n 

*nO,*nl,*n2,“‘,*nn  ^[0 


mi! 

m2! 

*’io!tn!ti2! . 

.MlJ 

*20^21  ^22!  •  • 

M2n\ 

m3! 

mn\ 

*3o!*31  !t32!  • 

^nO  *^nl  *^n2  •  •  •  • 

linn'- 

_*10_*U  «,*12  , 
X10  X11  J'12 

«*ln 

xln 

~,*20  ^*21  ^,*22 
x20  X21  x22 

r*2n 

*  '  *  x2n 

~,*30  *31  ^v, *32 

x30  X31  x32 

~*3n 

x3n 

„1ti0t*r1  rln2  . 

’  *  xn0  xnl  xn2 

.  ,  ™*nn 

^  nn 

to!ti 

!t2!  •  •  -!t„! 

(io  +  ^1  +  ^2  -f  •  •  *  +  hi 

4-  n)! 
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where 


U  —  *1J  +  *2/  +  *3/  +  •  "  +  ini,  l  —  0, 1, 2, . . . ,  n. 


Examples  of  Simplex  Integrations  on  Simplex 

A  general  n  dimensional  simplex  has  n+1  ordered  vertices  PqP\P2  . . .  Pn,  its  Jaco¬ 
bian  is 


X10 

X20 

X30  • 

. .  xno 

Xll 

X21 

Z31  . 

•  •  %nl 

X12 

X22 

X32  • 

*  •  %n2 

(15) 

n 

x2  n 

%Zn 

•  *  ^nn 

Formula  (14)  directly  gives  many  integrations  on  a  simplex. 


Sp0Pi  P2-Pn  (0, . . . ,  0,  mi,  0, . . . ,  0) 


HI  1. 


P 0  P 1  P2  •'  -  Pn 
-il2+ +  Ur 

**101*11 


x™1  du\du2duz  . . .  dun 


j  -  •  -piin  —  rni 

t  J  im.ii 


mil 


2 

X10  xl\  xl2 


X 


20!?i  !z2 !  •  ■  .HJ 


In 


(io  H~  i\  +  12  +  •  •  •  +  in  +  ft)! 

*io+*u “H/2+ — +*1 »» =mi 


T  v — ^iOT*iiT*f2T«"“r*in=w*i  •  .  .  •  rail 

j\  1 

^ j  im.il  1  .2 


xio  xli  x  12 


X 


In 


(mi  +  n)\ 


As  special  case,  the  normal  one  dimensional  integral  formula  (16)  can  be  derived. 
Sp0pl(m)  =  sign(J)  [  x™dxi—j'S^ 

J Pn  Pi  ^UO,»ll>C 


l  —m  rii 1 

Jiojll  _ 


"10  ^11 


(m  +  1)! 


1 


(»n  -  ®io) 


(m  +  1) 

(^n^io  +  xTi  1  x\o  +  xn  2 x\o  P  ■  ■  ■  P  x\ixTo  1+xiixio) 


m  + 


ilAll  ^lO  I 


(16) 


5PoPiPa(M)  =  sign(J)  J  J 


dx\dx2 


P0P1P2 
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SPoPiPaCM) 

‘S'PoPlP2(1>0) 

SPoPiPaUi0) 

Sp0Pi  P-2  (0»  1) 

>S'PoPlP2(2)0) 

‘5,p0PiP2(2»°) 

Sp0p1p2{Q,2) 


-'JL 


du\du<i 

UqUxU? 

uo+*n+U2=o  0! 


Ewii.^  (0  "I"  2)! 

=  sign(J)  /  /  dx\dxi 
J  J  Pn  P\  P'2 


=  \J 


sign(J)  /  / 

J  JP0P1P2 


EU0+U1+U2— 1  l! 

•  *  •  \n  TTXoYi^10  ^12 

Uo,*ii)*i2>0  (I  +  Z): 


^*10^*11  ^*12 
1  *^1 


=  sign(J)  /  / 

V  ^PoPx 


x\dx\dx2 


Pi 


(17) 


=  -J(x10  +  +  X12) 

=  sign(J)  /  /  x2dx\dx2 

J  J  Po  Pi  Pa 

=  ~J(X20  +  X21  +  £22) 

0 

=  sign(J)  /  /  x]dxidx2 

J  JPoPlPl 

y>*10+tU'f>12~2  2‘  TnorHlr»12 

=  J^uo,i„,M2>o  (2 +  2)! 3:103:11X12 

=  sign(J)  I  /  x\dx\dx2 

J  JPoPlPl 

=  —  J(Xj0  +  2T2!  +  X22  +  Xio^n  +  X10X12  +  ^11X12) 
1 2 


(18) 


J 


24 


(2xioXio  +X10X11  +2:10X12 

+X11X10  +2xiixn  +X11X12 
+X12X10  +X12X11  +2x12X12) 


=  sign(J)  / 

J  J  Pq  Pi  1 


x\dx\dx2 


=  — -  J(x20  +  X21  +  X22  +  X20X21  +  X20X22  +  X21X22) 

-“5 
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(2x20^20  +£2C>£21  +X2o£22 
+£21*20  +2x2i£2i  +£21*22 

+*22*20  +*22*21  +  2x22*22) 

Sp0p1p2(l,l)  =  sign(J)  /  x1x2dx1dx2 

J  Jp0p1p2 

\*10+*U+*12=1  v— ^i20  +  i21+i22=l 


J  ^^\*10-r*ll-t-2l2— -l  ^—vi20“ 

ho, in  ,*12>0 


lr*10rUl  r*12  Iri20_t2i  1 22 
•^10  X11  x12  1^20  ^21  x22 


'SPoPift(M)  =  <sz#n(J)  /  / 

</  ^P0PiP2 


>*21 1*22  >0 

(*10  +  ?2o)!(?11  +  *~2l)!(&‘l2  +  ?22)! 
(2  +  2)! 


—  2^^jXiqX2°  +  2*11*21  +  2x12X22 

+  *10*21  +  X10X22  +  £11X20  +  £11  £22  +  £l2£20  +  £l2£2l) 
=  + 


24 

(2xio£20  +£io£21  +£10X22 

+£ll  £20  +2xhX21  +£ll£22 

+£l2£20  +£l2£21  +2£i2£22) 


‘5,PoPiP2P3(0)0,0)  =  sign(J)  /  /  /  dxidx2dx3 

J  J  JP0P1P2P3 


(0  +  3)!  6' 

■SPoPiftPaUiM)  =  sign(J)  /  /  /  x1dx1dx2dx3 

J  J  JP0P1P2P3 

v*10+*n  +  *12  +  *13=l 


1*12 1*13^0 


-*10^*11  ^*12  *13 
^10  ^11  x12  X13 


1! 


'S,pop1p2p3(l,0,0)  =  sign(J)  /  /  /  x1dx1dx2dx3 

J  J  JP0P1P2P3 


(1  +  3)! 


—  ^ 24^Xl°  Xn  +  *12  +  £13) 


*S'PoPiP2P3(0, 1,0)  =  sign(J)  /  /  /  x2dxxdx2dx3 

J  J  J  Pq  Pi  P2  p3 


°^24^a:20  2:21  2:22  a:23) 


SpoPlP2P3(0,0,l)  =  sign(J)  /  /  x3dxidx2dx3 

J  J  ip0PiP2P3 


S’Pof’if'iPs^OjO)  —  sign(J)  ///, 


x\dx\dx2<lxz 

P0P1P2P?. 

_  j  y^!io  +  in  +  ti2+ii3=2  ^i12^ii3  . 

^“^*10  »*11  »*12>*13>0 


2! 


llOJll  112  r  *13  _ _ _ 

xio  xn  xi2  xi3  ^  +  3)1 


—  +  xh  +  xn  +  X13 

o(J 

+  ^11^12  +  X\1X13  +  ^12^13) 


+Xio£ll 


•S’PoPj  p2p3(2,  0)0)  —  sign(J) 

=  Jl20 

(2xio^io 
+£n£io  +2xxiXn 
+£i2£io  +Z12E11 

+£l3£l0  +^13^11 

SPoPiPzPs  (0>2,0)  =  sign(J) 


Po  Pi  P*  Pa 


:13  +  a^io^ll  +  x10x12  +  x10x13 

0 

x\dx\dx2dxz 


+X10X12 
-\-Xl\X\2 
+  2xi2^12 
+X13X12 


+  X10X13 
+X11X13 
+  X12X13 

+  2x13X13) 


P0P1P2P3 


x\dx\dx2dxz 


J- 


1 


120 

(2X20^20  d~x20x21  +X20^22 

+X21X20  +2X21X21  +X21X22 

-\-X22x20  ~\~x22x21  +2x22^22 

+X23X2Q  +  x23x21  +x23x22 


SpoP, P,P,(0, 0, 2)  -  sign(J)  J  J  J  ^ 


+  X20^23 

+ X21 X23 

+X22^23 

+  2X23X23) 


x\dx\dx2dx3 


=  J 


1 


120 

(2x30X30 


+X30X31 


+X31X30  +2x31X31 

+x’32a:30  +X32X31 

+X33X30  +  X33X31 

Sp0PiP2P3{li  M)  =  sign(J) 


+  X30X32 
+  ^31X32 

+  2X32X32 
+  0:33^32 


+  X30X33 
+  X31X33 
+  £32£33 
+  2X33X33) 


x\x-2dx\dx2dxz 

P0P1P2P3 

i4'*llH“*12“t"*13:::=  1  x  ’v*20~H*21  4"  *22  “I"  *23  1 

“^*20j*21  i*22j*23^0 


j  s*10"r*il"r*i2T*i3- 

^  *10 1*11 1*12 1*1 3  >0 

*10*11  *12  *13  r*20  *21  i  22  £*23 

£10  £11  £12  X13  ^20  X2\  X22  x23 

(iio  +  *2o)K*11  +  *2l)K*12  +  J22)Kil++i23)! 
(2  +  3)! 


Spop1p2p3(IA,0)  =  sign{J)  f  l  f  x1x2dx1dx2dx3 

J  J  J  Po  Pi  P2  P3 


(23) 
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=  J 


1 

120 


Spop1p2p3(1,0,  1) 


(2a;io3;20 

+£io£21 

+  X10X22 

+£io£23 

+X11X20 

+2xhx2i 

+£ll  £22 

+Xh£23 

+£  1 2  X20 

+  £12£21 

+  2X12X22 

+£12£23 

+3:13X20 

+  3213  £21 

+£13£22 

+2X13X23) 

=  sign(J) 

If  x\x^dx\dx2dx^ 

J  j  JPoPlPlPz 

=  J— 
120 

(2x10X30 

+£io£3i 

+£l0£32 

-f  £10^33 

+^11^30 

+2xnX3i 

+£ll£32 

+  ^11^33 

+£’12£30 

+£12£31 

+2xJ2£32 

+^12^33 

+£13X30 

+£13X31 

+£l3£32 

+2^13^33) 

=  sign(J) 

III  X2Xzdx\dx2dx$ 

J  J  d  Pq  P\  P2  P& 

=  J— 
120 

(2X20^30 

+£20£31 

+£20£32 

+£20£33 

+X21 X30 

+2X21X31 

+  £2l£32 

+X2l£33 

+X22^30 

+  £22£31 

+2x22£32 

+£22£33 

+£23^30 

+  £23£31 

+£23£32 

+2x23  £33 ) 

(24) 


One  dimensional  simplex  integration  (16)  is  same  as  conventional  integration.  Formulae 
(17)-(20)  are  2-d  simplex  integrations  over  triangles.  Formulae  (21)-(24)  are  3-d  simplex 
integrations  over  tetrahedrons. 


Simplex  Integration  on  General  Two-dimensional  Blocks 

Since  simplex  integrations  always  have  the  Jacobian  J  as  factor  and  J  is  an  oriented 
area,  the  integrations  on  positive  area  and  negative  area  can  be  neutralized.  Denote 


P1P2P3  '  '  '  Pn  Pn- 1-1  —  Pi 


with 


Pt  =  (xU,X2i) 


as  a  polygon,  the  integrations  (17)-(20)  on 


P1P2P3  -  •  •  Pn 
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are  computed  by  summations.  (25)-(28)  are  the  integrations  of  DDA  method.  Here  P0  can 
be  any  point,  Pq  =  (0, 0)  are  chosen  in  order  to  have  simpler  formulae. 


J  J '  dXldx2  =^2k=lSpopkpk+1(0,0) 


1  v — 


_  i  xi  k  x2  k 

2  Z—‘k= l  x\  jb+i  x2  i 

EX\ 

SpoPkpk+1(l,0) 


X1  fc+1  x2  fc+1 


f  [  11  1  Y''”  Xl  k  X 2  k 

I  /  xidx\ax2  =  —  > 

J  J(A)  6  ^k= 1  Xi  fc+i  X2  fc+l 

//  *,<&,<(*,  =  i  V”  11  *  I; 

J  J  (A)  6  ^k=  1  X!  fc+1  X2 

l  J{A)Xidxidx2  ~  X^fc  i^p*pw(2’C) 

_  J+  : 

19  2—-'fc=i  xi  u.1  X' 


®1  fc  2:2  fc 
fc+1  fc+1 


(#1  fc  +  ^1  fc+l) 


(^2  fc  *^2  fc  +  l) 


_  Y^n  X\  k  X2  fc 
“  12  ^fc=l  Xi  fc  +  l  X2  fc+l 

{X1  k  +  x\  fc+1  +  Xl  fcXi  fc+l) 


[  [  x\dx\dx2  =  — 

J  j(A)  1  1  2  24  ^*=1 


_  Jy  y^n  xi  ^  x2  fc 
24  Z^fc=i  fc+i  x2  fc+i 

(2xi  fcXi  fc  +xi  fcXi  fc.fi 

+xi  fc+ixi  fc  +2xi  fc+iXi  fc+i, 


f  f  x~dx\dx2  =  — 

J  J(A)  2  1  2  24  ^k=i 


_  _j_  fc  x2  fc 

24  ^fc=i  xi  fc+i  x2  fc+i 

(2x2  fc-x2  fc  +x2  fcx2  fc+i 
+x2  fc+1^2  fc  +2x2  fc+i x2  fc+i) 


J  j(A)XlX  2dXldX2=T,l=  ;1  *S'p0pfcpJC+1(l, !) 

_  J_  v-n  a-'i  fc  i 

24  ^k= 1  xi  fc+i  x- 


_  v-n  a,'i  fc  x2  fc 
~  24  ^fc=i  xi  fc+i  x2  fc+i 

(2xx  fcx2  fc  +  2xi  fc+i x2  fc+i  +  xi  fcX2  fc+i  +  xi  fc+i x2  fc) 


I  [  xix2dxidx2  =  7^7  Y\ 

J  JiA)  24  ^k= 1 


_  v-n  aii  fc  x2  fc 
24  ^fc=i  xi  fc+i  x2  fc+i 

(2xi  fc-x2  fc  +xi  fc-x2  fc+i 

+‘'ci  fc+i-T2  fc  +2xi  fc+ix2  fc+i) 
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Figure  3.  shows  equal-lateral  2-dimensional  polygons.  The  equal-lateral  triangle, 
square,  pentagon  and  hexagon  are  in  the  unit  circle.  The  distance  from  each  node  to  the 
circle  center  is  1.  The  edge  length  and  area  of  each  equal-lateral  polygon  are  listed  in  the 
following  table.  The  centers  of  gravity  of  the  equal-lateral  triangle,  square,  pentagon  and 
hexagon  are  the  center  of  the  unit  circle. 


Figure  3.  Equal-lateral  of  2-dimensional  polygons 


The  coordinates  (a;,-,  y.)  of  vertices  of  equal-lateral  polygon  can  be  computed  by  the 
following  formulae: 


Xi  =  sm(360/n  *  i) 
yi  =  cos(360/n  *  i) 
i  =  1,2,  •  •  •  ,  n 

here  n  is  the  edge  number  of  the  polygons. 

Using  the  simplex  integrations,  the  same  area  and  center  of  gravity  are  obtained  for 
each  equal-lateral  polygon. 
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polygon 


edge  length 


area 


equal-lateral  triangle 
equal-lateral  square 
equal-lateral  pentagon 
equal-lateral  hexagon 


x/3 

V2 _ 

1 


\3s/3  =  1.299038 
2.000000 

§V/10  +  2V/5  =  2.377641 
=  2.598076 


Simplex  integrations  on  General  Three-dimensional  Blocks 

Integrations  (21)-(24)  are  the  integrations  of  tetrahedrons.  By  summations  of  sim¬ 
plex  integrations  (21)-(24),  the  volume  or  integrations  of  any  3-d  block  can  be  computed. 
Assume 

pM  pH  p[*l . . .  pN 
M  ^3  rn(t) 

pH  _  pM  »  _  1  2  3  « 

Fn(i)+ 1-M 

are  all  outward  rotated  polygons  of  a  block, 

pbl  _  fT-bl  1 

rj  —  'L2j3  x2j! 

and  P0  =  (0, 0, 0).  The  volume  of  this  block  is  given  by  (29).  Computed  by  simplex  in¬ 
tegrations,  integrals  (29)-(32)  are  represented  by  the  coordinates  of  the  boundary  vertices 
only. 


[  f  [  dx\d,X2dxz  =  V"\  Sp  p[.]p(.]pW  (0,0 

J  J  J(V)  ^i=1  ^ k=1  PoPl  k  k+1 


,0) 


n(») 


6  ^ i  =  1 


'l  1 


'1  k 


X 


2  1 


y2  k 


^3  1 


X 


3  k 


X 


IN 

1  k- fl  *+1  *0  Jfc+1 


/ /  J(v)x'dxidx2dx3  i=i 

x 


l  l 


X 


X 


1  fc 
[*] 

1  fc-f- 1 


X 


"2  1 


"2 

«] 

2  Jt+1 


^3  1 


x3  fc 
[i] 

r’  L  J  . 

"3  ifc  +  1 


(29) 
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IIL 


(zi  1  +  Xy  k  +  Xi  k+i) 


X2dxidx2dx$  = 


1  \  '  "V5 

Z^i=i  Z-j k=\ 


24 

[«] 


4*1 


[*1 


X1  1 

T4l 

/yi  l  J 

X2  1 
f -i 

^3  1 
r  *i 

rW 

X1  fc 

x2  A: 

Til 

X3  fc 
Til 

,ilJ 

l  fc+l 

^LZJ 

^2  fc+1 

«»L2J 

x3  fc+1 

IIL 


(X2  1  +  X2  fc  +  ^2  fc+l) 

1  ^ ^ tt  (  Z  ) 

x3dx1dx2dxz  =  —  >  > 

24  'i=l  *-^k= 1 


24 

Tra 

xi  j 

xi  * 

j«] 

a  fc+i 


4*i 


[*] 


in 

m 


x\dx\dx2dx$ 

x\dx\dx2dx$ 


=  ES 

1 '*=i 


i 

120 

J*1 

xi  i 

J*1 

xi  jt 

w 


/*>■*  L  J 
x2  1 
r;l 

rp  l  J 

X3  1 

r -i 

x2  fc 

\i) 

rw 

x3  fc 
fil 

x2  fc+1 

^»L*J 

^3  fc+1 

fc  +  #3  fc+l) 

y"<0  s 

^fc=  1 

PoP^pJ 

■>  5  \  -\  n 

*^t=l  ^“—/fc 

(0 

=1 

J«] 

,c2  1 
r  -i 

x3  1 

r  *i 

X2  fc 
ni 

r  M 

X3  fc 

rvi 

<7*  L2J 

^2  fc+1 

x3  fc+1 

fc+1 


IIL 


x\dx\dx2dxz  =  - 


fc+i 

(2xi  ia?i  i  +xi  ixi  * 

+#i  fc^i  l  +2#i  k 

+#1  fc  + 1^1  1  +  ^1  fc+l^l  A; 

1  ^"■v5  v  ^n(i) 

=1  Z—/k=l 


+  Xj  \X\  Jfc+i 
+£l  k% i  fc+i 
+2xi  fc+ii'i  jk-t-i 


IIL 


x\dx\dx2dx$  = 


120  ^i=i 

J*1 

xi  i 

J*1 

xi  it 

4*1 

'1  k+ 1 
(2X2  1X2  1 
+X2  kX2  1 
+X2  fc+1^2  1 

1  V  ■  "V  3 
120 
J*1 


4*1 


[*] 


*21 
r -i 

/y>  1  J 

x3  1 

JlJ 

X2  fc 
r,i 

TllJ 

X3  fc 

fi’l 

,l?J 

2  fc+1 

xLlJ 

x3  fc+1 

+X2  1  X2  k 
+  2x2  k% 2  k 
+#2  fc+1^2  k 

2-^k= 1 


+  #2  1^2  fc+1 
+  ^2  fc#2  fc+1 
+  2^2  fc+ix2 fc+l) 


*1  1 
J*1 

X1  k 

[*] 

PL  J  , 

'1  fc+1 


X 


x2  1 

XM 

X2  k 

[*1 

2  fc  +  1 


x3  1 
r(*1 
X3  fc 

XH 

x3  fc+1 


(30) 
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(31) 


(2x3  1^3  1  +*3  1^3  i  +^3  1^3  HI 

+  X3  fc-X 3  1  +2X3  kX 3  k  +X3  kX 3  k+1 

+  X3  fc+iX3  i  +X3  fc+iX3  k  +2x3  t  +  lX3  k  +  l) 


UL 
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Figure  4.  shows  3  different  blocks.  Each  block  is  shown  from  two  different  view 
angles.  For  each  block,  the  geometric  formulae  of  boundary  plan  angle,  edge  length,  dis¬ 
tance  of  boundary  polygon  node  to  the  polygon  center,  distance  of  polygon  plane  to  the  unit 
sphere  center,  and  the  theoretical  volume  of  the  block  are  listed  in  the  following  tables. 

The  first  block  is  a  tetrahedron  with  four  equal-lateral  triangles  as  its  boundary  faces. 
The  block  is  in  the  unit  sphere.  The  distance  from  each  vertex  to  the  sphere  center  is  1. 
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First  block  of  4  face  polygons 


boundary  polygon  number 

edge  number  of  polygons 

angle  of  two  adjacent  polygon  planes 

edge  length  of  polygons 

distance  from  polygon  to  sphere  center 

distance  of  node  to  its  polygon  center 

distance  from  vertex  to  sphere  center 

volume  of  block 


4 

3 

arctan(  2\/2) 


l 

3 


fv/2 

1 

^  =  0.51320 


The  second  block  is  a  cube  with  six  equal-lateral  squares  as  its  boundary  faces.  The 
block  is  in  the  unit  sphere.  The  distance  from  each  vertex  to  the  sphere  center  is  1 . 


Second  block  of  6  face  polygons 


boundary  polygon  number 

edge  number  of  polygons 

angle  of  two  adjacent  polygon  planes 

edge  length  of  polygons 

distance  from  polygon  to  sphere  center 

distance  of  node  to  its  polygon  center 
distance  from  vertex  to  sphere  center 
volume  of  block 


=  1.53960 

3V3 


The  third  block  has  twelve  faces,  each  face  is  a  equal-lateral  pentagon.  The  block  is 
in  the  unit  sphere.  The  distance  from  each  vertex  to  the  sphere  center  is  1. 


Third  block  of  12  face  polygons 


boundary  polygon  number 

edge  number  of  polygons 

angle  of  two  adjacent  polygon  planes 

edge  length  of  polygons 

distance  from  polygon  to  sphere  center 

distance  of  node  to  its  polygon  center 
distance  from  vertex  to  sphere  center 
volume  of  block 


12 

5 

180  —  arctan{2) 

\ZE"+  2\/5 

TTs  V^2(5  ~  V/5) 

1 

|^y/10(3  +  \/5)  =  2.78516 


In  order  to  use  the  simplex  integration,  the  node  coordinates  (x,  y)  in  each  boundary 
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polygon  have  to  be  computed.  Based  on  the  angles  of  boundary  planes  and  distances  of  these 
planes  to  the  sphere  center,  the  polygons,  the  edges  and  the  vertices  can  be  computed. 

The  three  dimensional  simplex  integrations  give  the  same  volume  as  computed  di¬ 
rectly  from  the  geometric  formulae  listed  in  the  tables. 


Rock  Failure  Examples  by  DDA 

The  equilibrium  of  the  DDA  method  is  reached  by  minimizing  the  total  potential 
energy.  As  the  energy  is  computed  by  integrations,  most  of  the  DDA  formulae  are  formed 
by  the  polynomial  integrations  over  the  generally  shaped  blocks.  The  simplex  integration  is 
developed  and  applied  to  DDA  formulation.  This  new  integration  gives  analytical  solutions. 
The  integrand  can  be  multi-dimensional  polynomials. 

In  the  two  dimensional  case,  the  integration  domain  can  be  any  convex  or  concave 
polygons.  In  the  three  dimensional  case,  the  integration  domain  can  be  any  complex  body 
with  plane  polygon  boundary.  The  integration  results  are  simply  represented  by  the  coordi¬ 
nates  of  boundary  vertices.  Based  on  the  simplex  integration,  DDA  algorithms  are  simple, 
efficient  yet  accurate.  Most  important,  the  accurate  integral  solution  of  mass  matrix  ensured 
the  convergence  of  “open-close”  iterations. 

Three  rock  failure  examples  are  presented.  The  failure  process  is  a  transition  from 
continuous  to  discontinuous  states.  The  discontinuous  deformation  analysis  (DDA)  has  to 
fulfill  physical  laws  of  both  continuous  and  discontinuous  materials.  When  the  computed 
displacements  and  deformations  are  large  enough  to  be  visible,  the  mechanism  of  the  failure 
and  the  final  damage  can  be  shown,  and  the  ultimate  strength  of  materials  or  structures  can 
be  intuitively  estimated.  The  visible  sliding  and  joint  opening  from  the  computation  can 
demonstrate  that  the  physical  laws  are  satisfied. 

The  computations  require  equilibrium  in  both  the  discontinuous  contacts  and  the 
continuous  zones  throughout  the  entire  dynamic  process.  Following  a  real  time  sequence, 
the  DDA  uses  a  step  by  step  approach.  The  displacements  of  each  time  step  are  so  small  that 
normal  linear  equations  for  small  displacements  can  be  adopted.  At  the  end  of  each  time 
step,  the  equilibrium  in  both  discontinuous  interfaces  and  continuous  zones  are  reached. 

As  the  step  displacements  are  small,  the  kinematic  relation  and  friction  law  are  ex¬ 
pressed  as  linear  inequalities.  Based  on  natural  contact  phenomena,  an  “entrance  theory” 
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was  developed.  The  entrance  lines  are  used  to  form  linear  inequality  equations.  The  same 
linear  inequalities  are  used  to  define  the  entrance  distances,  entrance  points,  entrance  con¬ 
straints  and  entrance  criteria.  The  “open-close”  iterations  ensure  that  no  tension  and  no  pen¬ 
etration  occur  at  all  entrance  positions.  There  are  three  entrance  modes:  open,  sliding  and 
locking.  Coloumb’s  Law  is  also  fulfilled  at  all  entrance  modes  and  all  entrance  positions. 

There  are  1500  to  2000  rock  blocks  in  each  example.  The  dimensions  of  computed 
regions  are  about  40  to  80  meters.  The  numbers  of  time  steps  are  from  300  to  600.  The  total 
elapsed  times  are  from  0.2  to  2.0  second.  The  maximum  total  displacements  are  more  than 
ten  times  the  average  block  diameter. 

Figure  5  shows  the  collapse  process  of  a  tunnel  caused  by  high  initial  stresses. 

Figure  6  shows  the  penetration  of  a  missile  at  a  velocity  of  300  meters  per  second 
into  a  blocky  rock  mass  with  two  tunnels. 

Figure  7  shows  the  damage  state  as  a  strong  stress  wave  passing  through  three  tun¬ 
nels  excavated  in  a  blocky  rock  mass. 

Figure  8  shows  a  rock  toppling  failure  caused  by  slope  excavation. 
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Figure  5.  Collapse  of  a  tunnel  caused  by  high  initial  stresses  (continued) 
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Figure  5.  (Continued) 
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Figure  5.  (Concluded) 


The  penetration  ot  a  missile  at  a  velocity  or  juu  meters  pei 
second  into  a  blocky  rock  mass  with  two  tunnels  (continued) 
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Figure  6.  (Continued) 


Figure  7.  A  strong  stress  wave  passing  through  three  tunnels 
excavated  in  a  blocky  rock  mass  (continued) 
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Figure  7.  (Continued) 


Figure  7.  (Continued) 
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Figure  7.  (Concluded) 
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Figure  8.  (Continued) 
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Figure  8.  (Continued) 
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